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FLOWFIELD-DEPENDENT  MIXED  EXPLICIT-IMPLICIT(FDMEI)  ALGORITHM 
TOWARD  DIRECT  NUMERICAL  SIMULATION  IN  HIGH  SPEED  FLOWS 


T.  J.  Chung 

Department  of  Mechanical  &  Aerospace  Engineering 
The  University  of  Alabama  in  Huntsville 


SUMMARY 


This  report  covers  the  results  of  the  research  on  new  concepts  and  formulations  aimed 
toward  direct  numerical  simulation(DNS)  dealing  with  high  speed  flows.  The  research 
was  motivated  by  the  fact  that  it  is  desirable  to  develop  a  CFD  program  which  can  be 
applied  to  all  speed  regimes,  both  compressible  and  incompressible,  both  viscous  and 
inviscid,  and  both  laminar  and  turbulent  flows,  ideal  for  shock  wave  turbulent  boundary 
layer  interactions  in  hypersonics.  The  popular  notion  that  DNS  will  resolve  all  turbulent 
microscales  can  be  applied  to  incompressible  flows.  For  compressible  flows  with  shock 
waves  interacting  with  turbulent  boundary  layers,  however,  difficulties  arise  in  dealing 
with  complex  physical  phenomena  such  as  transition  from  laminar  to  turbulent  flows, 
relaminalization,  interactions  between  viscous  and  inviscid  flows,  and  high  temperature 
gradients  close  to  the  wall,  particularly  in  hypersonics.  No  currently  available  CFD 
techniques  are  capable  of  resolving  these  physical  phenomena  simultaneously  even  in 
DNS. 


The  purpose  of  the  FDMEI  approach  is  to  overcome  these  difficulties  by  introducing 
the  flowfield-dependent  implicitness  parameters  which  are  calculated  from  changes  of 
Mach  numbers,  Reynolds  numbers,  Peclet  numbers,  and  Damkohler  numbers(if  reacting) 
between  adjacent  nodal  points  and  between  time  steps.  These  implicitness  parameters 
imply  current  flowfields  changing  in  space  and  time  designed  to  alter  the  magnitudes  of 
every  term  in  the  Navier-Stokes  system  of  equations,  reflecting  the  parabolic,  hyperbolic, 
and  elliptic  nature  of  the  actual  flowfield.  For  example,  far  away  from  the  wall,  the  initial 
form  of  the  Navier-Stokes  system  of  equations  automatically  changes  into  a  hyperbolic 
form  of  the  Euler  equations  as  dictated  by  the  flowfiel-dependent  implicitness  parameters. 
By  the  same  token,  the  viscous  terms  become  activated  as  boundary  layers  are  approached 
and  they  become  dominant  close  to  the  wall,  again  automatically  dictated  by  the  flowfield- 
dependent  implicitness  parameters.  Such  phenomena  can  be  clearly  observed  when  the 
contours  of  these  fiowfield-dependet  implicitness  parameters  are  plotted,  which  are  shown 
to  be  representative  of  the  flowfield  themselves. 


The  essence  of  the  FDMEI  scheme  is  to  follow  the  physics  and  as  a  result  the 
numerics  are  generated  accordingly.  This  is  contrary  to  all  other  existing  computational 
schemes  in  which  numerics  are  predetermined  for  the  fixed  physics  under  investigation. 
Unfortunately,  howevser,  physics  change  significantly  as  a  function  of  space  and  time  for 
which  the  predetermined  numerics  are  no  longer  suitable.  This  occurs  when  the  program 
desinged  for  incompressible  flows  is  to  cope  with  compressible  flows  in  different  regions 
of  the  domain  and  vise  versa,  or  for  laminar  flows  to  handle  turbulent  flows  and  vise 
versa..  Such  computational  schemes  will  not  be  successful  even  in  DNS  mech 
refinements.  The  emphasis  of  FDMEI  is  upon  not  only  the  ability  to  deal  with  all  situations 
of  fluid  dynamical  physics  but  also  the  computational  accuracy  if  and  when  the  computer 
is  available  for  DNS  calculations.  In  the  mean  time,  the  FDMEI  approach  can  be  used  for 
non-DNS  problems  with  accuracy  much  superior  to  any  other  CFD  methods  available 
today. 


The  results  of  the  research  on  FDMEI  are  summarized  in  three  journal  publications 
attached  herein. 

1 .  Yoon,  K.  T.  and  Chung,  T.  J.,  “Three  Dimensional  Mixed  Explicit-Implicit  Generalized 
Galerkin  Spectral  Element  Methods  for  High-Speed  Turbulent  Compressible  Flows”, 
Computer  Methods  in  Applied  Mechanics  and  Engineering,  Vol.  135,  pp  343-367,  1996. 

2.  Chung,  T.  J.,  “A  New  Computational  Approach  with  Flowfield-Dependent  Implicitness 
Algorithm  for  Applications  to  Supersonic  Combustion”,  in  Avanced  Computational  and 
Analysis  of  Combustion,  Ed.  G.  D.  Roy,  S.  M.  Frolov,  and  P.  Givi,  Moscow;  ENAS 
Publishers,  pp.  466-489,  1997. 

3.  Yoon,  K.  T.,  Moon,  S.  Y.,  Garcia,  S.  A.,  Heard,  G.  W.,  and  Chung,  T.  J.,  “Flowfield- 
Dependent  Mixed  Explicit-Implicit(FDMEI)  Methods  for  High  and  Low  Speed  and 
Compressible  and  Incompressible  Flows”,  Computer  Methods  in  Applied  Mechanics  and 
Engineering,  Vol.  148,  1997. 

As  a  consequence  of  this  research,  the  following  important  conslusions  and 
recommendations  are  provided:  (1)  As  shown  in  Appendix  A  of  Reference  3  above,  the 
FEMEI  scheme  leads  to  all  existing  computational  methods  if  the  flowfield-dependent 
implicitness  parameters  are  fixed  to  certain  arbitrary  numbers  between  zero  and  one, 
indicating  that  all  existing  methods  are  the  special  cases  of  FDMEI,  (2)  The  FDMEI 
scheme  provides  a  single  computer  code  which  can  be  applied  to  all  physical  phenomena 
in  fluid  dynamics,  (3)  Due  to  limited  computational  resources  and  limited  number  of 
research  personnel,  benchmark  validations  included  only  small  portions  of  CFD  problems 
in  this  report.  The  future  research  should  include  higher  Mach  numbers  and  Reynolds 
numbers,  detailed  studies  of  transition  to  and  from  turbulent  flows,  high  temperature 
gradients,  compressibility  effects,  dilatational  thermal  dissipation,  and  finally  the  firm 
establishment  of  the  FDMEI  technology  benefiting  the  CFD  community  in  general. 
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Modeling  of  turbulence  has  been  the  controversial  subject.  Closure  models,  probability  density 
functions  (PDF),  large  eddy  simulation  (LES).  direct  numerical  simulation  (DNS)  and  other  methods 
have  been  reported. 

The  purpose  of  the  present  study  is  to  introduce  an  unstructured  adaptive  spectral  element  method  in 
dealing  with  combined  turbulence  and  shock  waves  for  both  internal  and  external  flows  of  aerospace 
vehicles.  This  work  is  motivated  by  the  fact  that  DNS  can  be  achieved  via  adaptive  h-p  methods, 
combining  the  mesh  refinement  (/i-method)  with  spearal  polynomial  degrees  of  freedom  (p-method)' 
It  is  well  known  that  the  most  crucial  aspect  of  turbulent  flows  is  microscales  involved  in  boundary 
layers  (viscous  sublayer,  buffer  zone  and  turbulent  core)..  This  is  where  the  spectral  polynomial  degree 
of  freedom  can  be  increased  as  desired  since  the  mesh  refinement  alone  is  incapable  of  resolving  the 
microscale  requirements..  In  this  way,  turbulence  modeling  techniques  can  be  avoided.  Furthermore, 
the  current  practice  in  DNS  to  use  extensive  refinements  in  finite  difference  discretization  may  also  be 
avoided.  Babuska  and  his  co-workers  [10-13)  and  Oden  and  his  co-workers  [14-17]  contributed  to  the 
advancement  of  FEM  h-p  adaptive  methods.  Their  applications  have  not  been  extended  to  shock  waves 
interacting  with  turbulent  boundary  layers.  In  what  follows,  the  p-version  of  the  finite  element  method 
ts  referred  to  as  the  spectral  element  method  with  Legendre  polynomials.  Although  the  terra  ‘spectral 
element  method'  was  used  by  other  investigators,  the  basic  approach  in  the  present  studv  is  sienificantlv 
different  from  the  earlier  work  [25-27  j. 

Chung  and  his  co-workers  [18-24)  have  studied  finite  element  strategies  as  applied  to  shock  wave 
turbulent  boundary  layer  interactions  m  non-reacting  and  reacting  flows.  The  main  emphasis  in  the 
present  study  is  to  establish  the  basic  theory  and  computational  strategies  of  MEI-GGM  involved  in  the 
Legendre  polynomial  spectral  element  method  and  to  present  preliminary  computational  results. 
Development  of  theory  and  formulation  include  irregular  node  connectivity  of  Legendre  polynomials  of 
various  orders.  Comparisons  with  experimental  results  have  demonstrated  superiority  of  the  direo 
numerical  simulation  over  the  standard  K-e  model  with  compressibility  effects  [19. 20).  One  of  the  most 
important  aspects  of  the  proposed  method  is  the  mixed  explicit-implicit  (MEI)  scheme  in  which 
flowfield  dependent  implicitness  parameters  as  calculated  from  local  Mach  number  and  Reynolds 
number  provide  automatic  adjustments  of  computational  processes  required  for  compressible  and 
incompressible  flows  or  high  speed  and  low  speed  flows.  ITtis  is  in  contrast  to  other  computational 
schemes  in  which  the  hyperbolic-elliptic  pressure  equation  must  be  solved  separately  to  provide 
pressure  correaions  when  the  compressible  flow  becomes  incompressible. 

In  what  follows  we  discuss  the  governing  equations  and  solutions  of  Navier-Stokes  system  of 
equations  via  Mixed  Explicit- Implicit  Generalized  Galerkin  Methods  (MEI-GGM)  in  Section  2  direct 
numerical  simulation  with  spectral  element  methods  in  Section  3,  calculations  of  DNS  perturbation 
variables  in  Section  4.  calculations  of  flowfield-dependent  implicitness  parameters  in  Section  5, 
numerical  applications  in  Section  6.  and  conclusions  in  Section  7. 


2.  Governing  equations  and  solutions  of  Navier-Stokes  system  of  equations 

A  convenient  form  of  governing  equations  for  compressible  viscous  flows  may  be  written  in  terms  of 
conservation  variables  as  follows: 

dU  aF,  dC, 

at  ax,  ^  dx,~^  (1) 

where 


~  p~ 

pv, 

0 

r  0  1 

u  = 

pv, 

pE 

pv,v,+pS,^ 
_{pE  +p)v,_ 

.  G,= 

- 1 

_ 1 

with  standard  definitions  given  by 
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The  solution  of  governing  equations  will  be  carried  out  using  the  generalized  Galerkin  approach  with 
test  and  trial  functions  given  by  isoparametric  and  Legendre  polynomials  by  means  of  mixed  explicit— 
implicit  schemes.  In  general,  explicit  schemes  are  inexpensive  but  less  accurate  in  comparison  with 
implicit  schemes  for  regions  of  high  pressure  or  velocity  gradients.  In  case  of  rapid  variations  of 
gradients  throughout  the  domain,  it  is  often  desirable  to  devise  a  scheme  in  which  implicitness  can  be 
adjustable  in  accordance  with  gradients,  more  implicit  for  the  region  of  high  gradients  and  less 
or  fully  explicit  for  the  region  of  low  gradients.  It  is  our  objective  to  obtain  amounts  of  explicitness  and 
implicitness  based  on  the  flowfield  dependent  parameters  such  as  TVD  limiters  in  FDM. 

To  this  end.  we  expand  in  Taylor  series  about  V'  by  introducing  the  implicitness  parameters  5, 
and  j,  for  the  first  and  second  derivatives  of  U  with  respect  to  time  (18-20],  respectively. 


t/-^' ={/'■  + A/  +2 

(2) 

with 

df/*”'  dU*  dSU"*' 

dt  ~  dt  ' dt 

(3a) 

■  at*  o«sj«i 

(3b) 

Substituting  (3)  into  (2)  yields 

fd^U”  a*A£/""'\ 

^  a<  j  +  o(A/>) 

(4) 

It  follows  from  (1)  that 

dU  _  dF,  dG, 

d/  ■  ax,  dx,  ® 

(5) 

Here,  F,  is  a  function  of  U  and  G,  is  a  function  of  U  and  its  gradient  t/*  so 
convective  Jacobian  c„  dissipative  Jacobian  b,  and  dissipative  gradient  Jacobian  c 

that  we  denote  the 

as 

_  BF,  BG,  BG, 

Note  that  if  the  source  teim  B  includes  variables  such  as  reaction  rates  then  it  would  be  necessary  to 
consider  the  source  term  Jacobian. 

The  second  derivative  of  U  with  respect  to  time  may  now  be  written  in  terms  of  these  Jacobians. 

Bt^  (“(  /)  dx,  (  Bxj  dx,  Bx,  dx*  (  BXj  dx, 

(6) 

Substituting  (5)  and  (6)  into  (4)  and  neglecting  the  third-order  spatial  derivatives  associated  with  c,* 
yield 

">  =  iy>.i  +  '^i.i  ~  J  ^k.k^il)  E  =  U.y, 
M  ^  T  +  Sq  \T  )  5o**110K 

p  =  pRT 


A£r"^'  =  M  - 
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dx,  ax, 


—  c  1 

(dAF"' 

>  5  J 1 

\  dx, 

'dF” 

oG" 

~~ 

{a.  +  b)  —  (—^  —  -RA  t  u.  ^  \1 

+  0(At')  '  ' 

In  order  to  provide  implicitness  to  diffusion  behavior  diffcrcnriv  t. _ 

5,  associated  with  G„  respectively,  as  follows  fferently  from  convection,  we  reassign  j,  and 

J,  AG,  ^  iiAG,  s^AG,  ^  s^AC, 

■  (8a,b) 

with  the  vanous  implicitness  parameters  defined  as 

5 1  —  first-order  convection  implicitness  parameter 

J;  =  second-order  convection  implicitness  parameter 

Sj  =  first-order  diffusion  implicitness  parameter  (8c) 

J4  =  second-order  diffusion  implicitness  parameter 

-£E 

,,  ,a=4t,- 


1 

a^AG"*'\ 

‘  ax,dx^.  ) 

•^2  2 

(a, a,  + 

A/'  r 

d 

/dF“ 

-— 1(». 

u 

'“*>  ta'  A2 

The  Galerkin  integral  of  (9)  may  now  be  carried  out  as  follows  (21]: 

/  <P„F(f/.F,.C,)dG=0 

(10) 

are  inteipolated  by  theS^fuSons'"0["i'^°"’*'^  functions,  and  the  conservation  variables 

U(x,  0  =  0Jx)UJO  F,(x,  /)  =  0Jx)F^XO  G,(x,  r)  =  0„(x)C,,(O 
Substituting  (II)  into  (9)  and  (10)  yields 


—  f 

2  ^ irqb dG 

+i{-K  <•.(*, -,„4  +c„„4  ,)| 
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HI  =  {st[0^_,cp,iF;.,  -f  g;,,)  + 

-  ^  + g;>j  -  0...0,Bi]  -  ^  0^0,Bi]  dn 

Nl  =  -  G" )  -  ^  (a,,,  +  6,„)(F;.,  ^  G;,;, dF 


« 

where  0^  indicates  the  Neumann  boundary  interpolation  function  with  unity  if  applied  and  zero 
otherwise.  Note  also  that  boundary  terms  in  are  assembled  only  into  boundary  nodes  of  B  .  For 
three  dimensions,  i,  j,k  =  1.2,3  associated  with  the  Jacobian  imply  directional  identification”©?  each 
Jacobian  matrix  (a,,  a,,  a,,  b,.  b,,  b,.  c,,.  c,,.  c,3,  c^,.  c,,.  c,3,  C3,,  £3,,  Cj,)  with  r.  j  =  1,2,3,4,5 
denoting  entries  of  each  of  the  5x5  Jacobian  matrices.  It  should  be  warned,  however,  that 
Jacobian  matrices  must  be  multiplied  precisely  as  dictated  by  summing  through  repeated  indices,  not 
through  matrix  multiplications  as  a  whole.  Indices  can  be  reduced  similarly  for  2*D. 

It  is  interesting  to  note  that  all  implicitness  parameters  can  be  shown  to  be  funnions  of  flowfields 
between  upstream  and  downstream,  and  that  the  covection  implicitness  parameters  r,  ands, 
with  the  first  term  in  are  analogous  to  the  total  variation  diminishing  (TVD)  limiters  in  the  FDM 
literature  isee  Appendix  A).  With  an  adequate  choice  of  these  implicitness  parameters,  acceptable 
resolutions  of  shock  waves  have  been  verified. 


On  the  other  hand,  the  diffusion  implicitness  parameters  s,  and  s/ are  capable  of  alleviating  and 
accommodating  the  stiffness  involved  in  turbulent  diffusion  or  finite  rate  chemistry  (if  reacting).  No 
analogy  can  be  shown  since  they  do  not  exist  in  other  numerical  schemes.  It  should  also  be  noted  that 
interactions  between  convection  and  diffusion  are  achieved  by  means  of  the  terms  associated  with  the 
products  and  These  terms  are  particularly  important  for  shock  wave  turbulent  boundary 

layer  interactions  where  the  effect  of  convection  upon  diffusion  and  vice  versa  is  crucial  in  order  to 
resolve  turbulence  microscales  as  disturbed  by  shock  wave  interactions.  We  shall  refer  to  these  terms  as 
convection-diffusion  interaction  terms. 

If  the  high  speed  compressible  flow  far  from  the  wall  becomes  the  low  speed  incompressible  flow  in 
the  vicinity  of  the  wall,  we  question  how  pressure  can  be  calculated  where  the  perfect  gas  law  is  no 
longer  valid.  To  this  end  we  integrate  the  steady  state  incompressible  momentum  equation 

\{p^\  pv,v^  ,  dx,  ^  f  (mu,,,,  +  )  dx, 


or 

1 

where  Py  is  the  constant  of  integration  and 

2  ~ ^  J  (ft  —  spatial  dimension) 


(13) 


where  lu,  is  the  component  of  vorticity  vector.  Eq.  (13)  as  related  to  the  perfect  gas  law  leads  to 


^0  =  pyCpT-  E  +  v,Vi)  -  Q 


(14a) 


or 


n  1 

Po  =  P«r  +  2  pw,u, -Q  .  ( 14b) 

If  the  constant  of  integration  or  stagnation  pressure  as  given  by  (14)  indeed  remains  constant,  then  this 
implies  that  the  compressible  flow  has  turned  into  an  incompressible  flow.  If  P-  as  calculated  from  (14) 
does  not  remain  a  constant,  then  the  incompressible  flow  has  been  changed  back  to  compressible  flow. 
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This  approach  allows  the  use  of  the  conservation  form  of  the  Navier-Stokes  system  of  equations 
througnout  the  domain  with  compressible  and  incompressible  elements  being  treated  accordinely  via 
flowfield  depenaent  implicitness  parameters  without  resonine  to  the  separate  hvperbolic-elliotic 
pressure  equation  for  pressure  correction  when  the  flow  becomes  incompressible.  ' 


3.  DNS  via  unstructured  spectral  element  methods  (h-p  version) 

Our  objective  here  is  to  resolve  time  and  length  scales  involved  in  turbulence  interacting  with  shock 
waves  using  adaptive  unstructured  h-p  finite  elements,  refered  to  as  spearal  element  methods  One 
approach  is  to  refine  the  mesh  (/i-methods)  until  further  refinement  is  unproductive,  at  which  time  the 
spectral  degrees  of  freedom  (p-methods)  are  increased  in  order  to  reduce  errors  as  desired  such  as  in 
the  region  of  turbulent  viscous  sublayer.  However,  the  more  desirable  approach  is  to  optimise  between 
the  mesh  refinement  and  spearal  orders.  Thus,  the  most  crucial  aspea  of  the  h-p  methods  is  to 

.  u"*  the  local  error  to  a  minimum 

Should  h  (mesh  size)  be  deaeased  or  should  p  (polynomial  or  spearal  degrees  of  freedom)  be 

maeased?  Although  some  work  in  optimization  between  the  h-  and  p-processes  has  been  reported 
1^4-161.  the  suDiect  of  opnmization  appears  to  be  an  open  question.  Thus,  our  approach  in  this  studv  is 
o  refine  the  mesh  until  shock  waves  are  adequately  computed  and  then  resort  to  the  p-version  with 
higher-order  Legendre  polynomials  in  order  to  resolve  turbulence  miaoscales.  Toward  this  end  the 

e^or  indicator  i)  may  be  defined  in  terms  of  density  for  shock  waves  and  velocity  for  turbulence.  Some 
of  the  options  are  given  as  follows: 


max 

s  =  hn^\p\Hil\p\„i 

max 

where  h  is  the  mesh  parameter  and  various  Sobolev  space  (H")  seminorms  are  defined  as 


dp 

dX: 


dn 


Ml 


IpI 


//j 


=  ([  ^'p  V'" 

Ufl,  dx,  dXi  dx,  dx^ 


d^V: 


d^u. 


dXj  dx^  aor,  dx* 


\W2 

dnj 


(15a) 

(15b) 

(16a) 

(16b) 


^e  choice  among  these  options  depends  on  various  physical  aspects  of  the  given  problem  whether 
dominated  by  density,  velocity  components,  their  gradients,  or  second  derivatives  For 

^sociated  with  shock  waves  and  the  p-adaptivity  associated  with  turbulent  boundary  layers,  respective- 

Direct  numerical  simulations  for  turbulent  flows  are  achieved  by  higher  snectral  orders  usine 

““'Polaiion’iuJion,  L 

TSg  1  S-o-te  isoparametric  element  are 


8  corner  node  isoparametric  functions: 
1 


‘*'v’  =  8(l+fvf)(l+77^T,)(l  +  ^^n 
12  edge  mode  Legendre  polynomial  functions: 

with  m  =  2 - ,q\  \2{q  -  1)  edge  modes;  q^2 


(17) 


(18) 
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6  face  mode  Legendre  polynomial  functions: 


(19> 


^ . 9-2:m■^„-4 . 9;3(,-2)(,-3)face  modes;  ,  »4 

1  interior  mode  functions  (bubble  functions): 

^''^lr=GM)G„(r,)G^U) 

with  m,n,p  =  2  /j-dm-i _ _  ^ 

„  .  . "  4,  .  .  „  =  6 . -  3)(,  -  4)(,  -  5)/6  intenor  modes;  ,  .6 

interpolation 


V2(2m-  i) 


(21) 

Similar  results  are  obtained  for  the  n  and  ^  H' 

dimensional  case  of  various  mixed  Legendre  oo[vmjm*iai  of  illustration  the  two- 

turbulent  microscale  distributions  is  shown  in  Fig^2^In  additio?^n?h  depiaing  a  possible 

asjace  1  (SI)  we  may  use  another  option  of Ve  ®  Polynomial  space  (defined 

m^es  and  (,  -  1)^  imerior  modes  (i,  1 2)  aJe  applIedTlif 

functions,  any'’ variable  U  may  b^ 


mmp^mmp 


where  U  ‘ 

integrals  in  addition  tiTTll);  ^  coefficients  to  be  determined  by  solving  the  following  Galerkin 


0lf‘/?(A£/)  d/2  =  0 

l^<t>Ly(Ac/)dn  =  o 
/^0L‘i^(A£/)d/2  =  O 


(23a) 

(23b) 


(23c) 
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q=l 

qmS 

qas? 

q=3 

qa2 

q=i 

qsl 

Rg.  2.  2-D 
modes). 


Interpolation  funaions  constniaed  by  Legendre  polynomial. 


<p<ci 

^  .V 


(comer  nodes).  (side  modes).  (intenor 


KnS„ 


^mkB^ri  ■*"  ^mkBn 

L^mkua^r,  +  ^mkum^r,  +  ^mkmnr. 


<nAs 


fl” 

anprs 


- 1 

w-*-! 

1 — 

1 _ 

Af/J. 

Ai/;. 

AU 

=  1^12; 

1— *6;  n 

AZ.s„*cz,„ 

^mkmp^rs  ^l^knprs 
^mkmmp^rs  ^Zkttmprg 


^artpq^rs  ^anpqrs 

^4. 

^mnpq^rs  ^  ^  mnpqrs 


^mkmpq^rs 


^  mknpqrs 


^mkumpq^n  ^mkumpqrs^ 


(24) 


modes:  a  =  comer  nnrip  ^  ^  ui  irccuom  rrom  edge,  face  and  interior 

integrals  are  shown  in  Appendix  B.  ’  ^  conservation  vanable  degrees  of  freedom.  Explicit  forms  of 

We  may  initially  consider  only  the  comer  node  equations. 
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+  cl* 


mknrs 

AlkuAs  +  ^mkunr. 


mnprs 

D^\ 

mknprs 


1 

3 

n 

= 

- 

Xikr 

^^mkur^ 

where 


^L,  =  (>lU5„  +  5l*,„)A£/3, 

^mkur  ~  ■*■  ^mku6rs) 


5  +C'' 

*  *  mnpq  rs  mnpqrs 

A  mknpq  ^rs  ^  mknpqrs 


A  q  X  4- 
^mkunp^rs  ^  mkunprs 


*■  mkunpq  ^rz 


5  ■+£ 


mkumpqrs^  J 


AC„ 


(26) 


which  act  as  source  terms  or  coupling  effect  of  the  comer  nodes  upon  spectral  behavior  through  edge, 
face  and  interior  modes.  The  final  step  is  to  combine  (25)  and  (26)  by 

At/;:'  =  -  i'a,  (27) 

with 


J'a,  =  (>4L5„  +  Blj  AUl  +  (/i:,,5„  +  J  At/:,.  -  (A^.pqS„  +  At/„„. 

Thus,  the  convergence  toward  shock  wave  turbulent  boundary  layer  interactions  can  be  achieved 
through  iterations  between  (26)  and  (27).  Note  that  in  this  process,  the  convection  implicitness 
parameters  s,  and  s,  are  held  constant  whereas  the  diffusion  implicitness  parameters  Sj  and  are 
updated  through  Reynolds  numbers. 

Our  objective  here  is  to  satisfactorily  simulate  turbulent  microscales  within  an  element.  All  edge,  face 
and  interior  mode  interpolation  functions  vanish  at  the  comer  nodes  but  exhibit  high  frequency 
variations  according  to  the  order  of  Legendre  polynomials  along  the  edges,  faces  and  interior  domain. 
It  is  intended  that  such  Legendre  polynomial  microscales  be  capable  of  simulating  the  physical 
microscales  of  turbulence  which  arc  involved  in  viscous  sublayer,  buffer  zone  and  turbulent  core.  The 
/i-adaptivity  alone  is  severely  limited  and  naturally  we  seek  a  remedy  of  this  situation  in  the  h-p 
adaptivity  utilizing  the  adequate  spectral  orders  required  for  accuracy.  Irregular  nodes  (hanging  nodes) 
which  arise  in  the  process  of  fi-adaptivity  are  treated  similarly  as  in  [14).  Furthermore,  the  advantage  of 
Legendre  polynomials  is  an  case  in  dealing  with  edge,  face  and  interior  modes  which  do  not  require 
specification  of  nodes  physically  located  in  the  element.  This  is  especially  beneficial  for  edge  and  face 
modes  in  establishing  boundary  continuities.  Continuity  of  variables  and  gradients  along  the  inter¬ 
element  boundaries  is  to  be  dictated  by  the  higher-order  polynomials  between  the  two  adjacent 
elements. 

For  two-dimensional  applications,  edge  and  face  modes  are  merged  to  side  modes  as  shown  in  Fig.  2. 
Consequently,  the  matrix  equation  (24)  can  be  reduced  so  that  only  side  and  interior  modes  arc 
retained. 


4.  DNS  perturbation  variables 

It  is  well  known  that  DNS  is  expected  to  provide  information  in  turbulence  microscale  levels  at  the 
expense  of  excessive  refinements  of  domain  discretization  [6].  The  purpose  of  the  present  study  is, 
instead,  to  avoid  such  refinements  by  means  of  implementing  high  spectral  Legendre  polynomial  orders. 
The  Navier-Stokcs  solver  as  introduced  here  allows  unsteady  time  accurate  solutions  from  which 
perturbation  variables_(/')  can  be  calculated  as  the  difference  between  the  Navier-Stokes  solution  (/) 
and  its  lime  average  /  [22,  23  j. 
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r-f-f  (28) 

This  computation  can  be  conducted  throughout  the  Navier-Stokes  integration  time  steps  or  upon 
arrival  at  quasi-steady  state.  Strialy  speaking,  in  shock  wave  turbulent  boundary  layer  interactions  a 
complete  steady  state  is  never  realized  as  unsteady  eddy  motions  persist  indefinitely,  althou^ 
background  flownelds  may  become  steady.  This  is  referred  to  as  the  quasi-steady  state.  The  time 
average  of  Navier-Stokes  solution  is  performed  using  the  Gaussian  quadrature.  In  this  process 
complicated  physical  phenomena  such  as  homogeneous  and  inhomogeneous,  isotropic  and  anisotropic, 
and  non-stationary  nature  of  perturbation  flowfields  in  shock  wave  turbulent  boundary  layer  interac¬ 
tions  can  be  resolved. 

Furthermore,  all  perturbation  variables  as  calculated  from  (28)  can  be  transformed  via  fast  Fourier 
transform  to  generate  power  spectral  density  vs.  frequency  domain.  Various  penurbation  variables  as 
well  as  background  flowfield  data  have  been  examined  [24].  As  a  result  of  this  study,  more  details  of 
shock  wave  turbulent  boundary  layer  interactions  such  as  variations  of  turbulent  kinetic  energy  vs. 
shock  strength,  laminar-turbulence  transition  instability,  relaminarization.  effects  of  dilatation,  etc.,  can 
be  rigorously  examined  in  comparison  with  the  previous  investigations  [6,  28.  29].  Some  limited  results 
and  discussion  of  these  subjects  are  presented  in  [24]. 


5.  Calculations  of  flowfieid-dependent  implicitness  parameters 


The  success  of  the  spectral  element  method  {h-p  version)  described  above  depends  on  accurate 
calculations  of  flowfieid-dependent  implicitness  parameters.  The  first-order  convection  and  diffusion 
implicitness  parameters  are  calculated  from  the  local  Mach  number  and  Reynolds  number  as  follows: 


s 


1 


’min(r.  1)  r^a 

-0  r<a,  Af„,„7^0  =  max(l  - s,, 0.5) 

.1 


with 


(29) 


r 


AAf 


(30) 


where  AM  is  the  difference  between  the  maximum  and  minimum  Mach  number  (AM  =  M„,„  - 
within  a  finite  element,  and  a  is  a  user-specified  small  number  (a  =0.01). 


h 


'min(s.  1)  J  3=  /3 
.  0  s</3,  Re„„7^0 


with 


5  = 


ARe 

Rc„.. 


Sf  -  max(l  -Sj,0.5) 


(31) 


(32) 


where  ARe  is  the  difference  between  the  maximum  and  minimum  Reynolds  number  (ARe  =  Re„„  - 
R®min)  within  a  finite  element  and  ^  is  a  user-specified  small  number  {p-0.01). 

The  flowfield  dependent  implicitness  parameters  as  defined  above  are  capable  of  allowing  various 
numerical  schemes  to  be  automatically  generated,  as  summarized  below: 

(1)  The  first-order  implicitness  parameters  s^  and  Sj  control  all  high  gradient  phenomena  such  as 
shock  waves  and  turbulence.  These  parameters  as  calculated  from  the  changes  of  local  Mach 
numbers  and  Reynolds  numbers  within  each  element  are  indicative  of  actual  flowfields.  The 
contours  of  these  parameters  closely  resemble  the  flowfields  themselves,  with  both  s,  and  s, 
being  large  (close  to  unity)  in  which  high  gradients  of  variables  exist,  but  small  (close  to  zero) 
where  such  gradients  are  small. 
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(2)  The  second-order  implicitness  parameters  s,  and  s,  are  also  flowfield  dependent.  However,  their 

primary  role  is  to  provide  adequate  computational  stability  (artificial  viscosity)  as  they  were 
originally  introduced  into  the  second-order  time  derivative  term  of  the  Taylor  series  expansion  of 
the  conservation  flow  variables  t/"”'.  Thus,  their  flowfield  dependency  is  limited  by  2=0,5 

for  adequate  computational  stability. 

(3)  The  ij  terms  represent  conveaion.  This  implies  that  if  5,  =0  then  the  effect  of  conveaion  is 
small.  The  computational  scheme  is  automatically  altered  to  take  this  effea  into  account,  with 
the  governing  equations  being  predominantly  parabolic. 

(4)  The  Si  terms  are  associated  with  diffusion.  Thus,  with  s,  =  0,  the  effect  of  viscosity  is  small  and 
the  computational  scheme  is  automatically  switched  to  that  of  Euler  equations  where  the 
governing  equations  are  predominantly  hyperbolic. 

(5)  If  the  first-order  implicitness  parameters  s,  and  53  are  non-zero,  this  indicates  a  typical  situation 
for  hyperbolic,  parabolic  and  elliptic  nature  of  the  Navier— Stokes  system  of  equations  with 
convection  and  diffusion  being  equally  important.  This  is  the  case  of  incompressible  flows  in  low 
speed.  The  unique  property  of  the  MEI-GGM  is  its  capability  to  control  pressure  oscillations 
adequately  without  resorting  to  the  separate  hyperbolic-elliptic  pressure  equation  for  pressure 
corrections.  The  capability  of  MEI-GGM  to  handle  incompressible  flows  is  achieved  by  a  Hpiirati. 
balance  between  s,  and  s,  as  determined  by  the  local  Mach  numbers  and  Revnoids  numbers.  If 
the  flow  is  completely  incompressible  (M  =  0),  the  criteria  given  by  (30)  leads  to  j,  =  1,  whereas 
the  implicitness  parameter  Sj  is  to  be  determined  according  to  the  criteria  given  in  (32). 


6.  Application 

As  benchmark  problems  two-dimensional  shock  wave  turbulent  boundary  layer  interactions  on  a  flat 
plate  and  compressible  comer  were  solved  in  [31]  and  demonstrated  an  excellent  comparison  with 
experimental  results  and  other  numerical  methods  (32],  Other  benchmark  problems  including:  (1)  the 
flat  plate  supersonic  boundary  layer  flow,  (2)  shock  wave  turbulent  boundary  layer  interactions  on  a 
compression  comer,  and  (3)  the  three-dimensional  sharp-leading  edged  fln  for  swept-shock  wave 
turbulent  boundary  layer  interactions  have  been  investigated,  which  are  presented  below. 

6.1.  Flat  plate  supersonic  boundary  layer  flow 

Rg.  3  shows  a  spectral  element  mesh  (gray  elements)  and  the  corresponding  density  contours  for 
Carter’s  flat  plate  problem  (33).  The  spectral  elements  appear  in  the  boundary  layer.  TTie  computed 
w^  pressure  and  skin  Motion  distributions  with  and  without  spectral  element  meshes  are  compared 
with  Carter  s  numerical  data  [33]  in  Fig.  4.  Note  that  symbols  SI  and  S2  imply  space  1  and  space  2 


I’ A****®^  eiemenu)  and  density  contoun  for  Carter’s  flat  plate  problem  with  boundary  conditions.  M,  -  3, 

T,  »  390  R.  =  2.8T..  Re^  -  1000. 
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0  02  Q4  QuB  QJ 


X 

(a)  Wall  pressure  profile  for  Imearai^  spectre!  elemem  nreshes 
conpannf  with  Gaiter's  data  (33]. 


X 


(b)  Skin  fhcuon  profile  for  linear  and  spectral  element  meshes 
coopanng  with  Cancr's  data. 

Fig.  4.  Companson  of  wail  pressure  and  skin  friaion  distribution  for  linear  and  spectral  meshes  with  Carter's  numerical  data. 


functions  of  Legendre  polynomials,  respectively,  as  defined  in  Section  4.  It  is  seen  that  for  the  laminar 
flow  the  spectral  element  method  does  provide  the  results  in  agreement  with  other  computational 
methods.  However,  the  spectral  elements  are  more  effective  when  the  flowficid  contains  large  gradients 
m  which  other  computational  methods  are  incapable  of  resolving  large  gradients  such  as  in  high 
Reynolds  number  and  high  Mach  number  turbulent  flows. 

6.2.  Two-dimensional  shock  wave  boundary  layer  interaction  on  a  compression  comer 

In  high  speed  vehicles,  deflection  of  a  control  surface  such  as  body  flaps,  elevons  and  rudders  causes 
the  interaction  of  shock  wave  with  boundary  layer,  which  may  cause  flow  separation,  resulting  in  a 
significant  decrease  in  flight  performance  and  excessive  increase  in  heating  rate.  A  two-dimensional 
compression  comer  experiment  by  Settles  et  al.  [35]  is  modeled  here  as  a  benchmark  problem  in  shock 
wave  turbulent  boundary  layer  interaction. 

Computational  geometry  and  scales  corresponding  to  the  experiment  are  shown  in  Fig.  5.  The 

stagnation  pressure  of  6.8  atm.  stagnation  temperature 
of  268  K.  freestream  unit  Reynolds  number  of  7.3  x  10 Vm,  deflection  angle  of  16«,  and  the  incoming 
turbulent  boundary  layer  thickness  of  2.3  cm. 

Adaptive  spearal  element  mesh  configurations  are  plotted  in  Fig.  6.  The  mesh  refinement  is 
performed  along  shock  waves  while  the  spectral  degree  is  increased  in  boundary  layer.  Convergence  of 
wall  pressure  for  different  Legendre  polynomial  spaces  and  degrees  at  :c/5  =  0.14  for  a  typical  transient 
state  IS  plotted  in  Fig.  7.  It  appears  that  the  convergence  rate  of  Legendre  polynomial  space  2  is  much 
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Fig.  5.  Compuiaiional  geometry  and  scales  for  16®  compression  comer.  The  freestream  conditions  corresponding  to  the  Settles' 
expenment  arc  M  -  2.85.  P,  -  6.8  atm.  T,  =  268  K,  Re,  =  7.3  x  lO’/m,  6^  =  2.3  cm. 


(b)  Focus  on  spectral  mesh  (gray  color)  in  boundary  layer. 


Fig.  6.  Adaptive  spectral  element  mesh  configuration.  The  mesh  is  refined  along  shock  waves  while  the  spectral  degree  is 
increased  in  boundarv  layer. 


Lagendre  polynomial  degree 

Fig.  7.  Convergence  of  wall  pressure  for  different  Legendre  polynomial  spaces  and  degrees  of  x/S  =0.14  for  a  typical  transient 
state.  It  appears  mat  the  convergence  rate  of  polynomial  space  2  (S2)  is  much  more  rapid  than  that  of  polynomial  space  I  (SI). 
Also,  the  higher  polynomial  degrees  are  closer  to  the  exact  solution,  independent  of  the  polynomial  spaces. 


(b)  Mach  number  contoura  (mtnaO.  maxs3  J8) 


(b)  contonis 

Fig.  8.  Density  and  Mach  number  contours  for  shock  wave  boundary  layer  interaction  on  a  compression  comer. 

■'"PUatness  parameters  (j,  and  s, )  contours.  The  s,  and  s,  contours  show  the  trend  of 


more  rapid  than  that  of  polynomial  space  1.  Also,  the  higher  polynomial  degrees  are  closer  to  the  exact 
solution,  independent  of  the  polynomial  spaces. 

Density  and  Mach  number  contours  are  shown  in  Fig.  8.  To  show  the  effect  of  implicitness 
parameters  upon  the  flowfield  calculations  the  contours  of  the  first-order  implicitness  parameters  s,  and 
are  plotted  m  Fig.  9  Note  that  5,  =0  away  from  the  shock  waves  and  boundary  layers,  but  becomes 
unity  at  locations  of  high  gradients  (shock  waves  and  boundary  layers).  By  the  same  token,  the 
impliatness  parameter  s,  representing  diffusion  behaves  similarly,  being  zero  and  unity  at  locations  of 


Fig.  10.  Convergence  history  of  energy  vanable  for  the  MEI  and  explicit  schemes, 
much  more  rapid  than  the  explicit  scheme. 


The  convergence  rate  of  the  MEI  scheme  is 
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Fig.  l..  Compamon  of  mean  streainw.se  velocity  profiles  for  the  present  results  and  Settles'  expenmental  data  at  several 
str^.nw,se  stauo.«.  The  plots  tllustrate  the  changes  ,n  boundary  layer  veloctty  profiles.  The  figure  at  x/6,  =  1.6  is  the  incoming 
cqmhbnum  turoulent  bounda^  layer.  The  location  of  the  shock  wave  appear  as  a  kink  is  some  velocity  profiles  downstream  of 

the  compression  comer.  The  downstream  profiles  are  seen  to  recover  rapidly  from  the  retarding  effects  ot  the  imposed  adverse 
pressure  gradjcnis.  ^  ^ 
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low  gradients  and  high  gradients,  respectively.  Note  that  the  contours  of  first-order  imolicitness 
parameters  j,  and  s,  resemble  the  flowfield  itself.  Here,  the  second-order  implicitness  parameten  5 
and  althougn  flowfield  dependent  to  some  extent,  their  pnmary  role  is  to  assure  computational 
stability  with  their  criteria  given  in  (30)  and  (32).  The  convergence  historv  of  energy  variable  for  MEI 
and  explicit  schemes  IS  compared  in  Fig.  10.  The  convergence  rate  of  the  MEI  scheme  is  much  more 
rapid  than  the  explicit  scheme.  Fig.  11  shows  a  comparison  of  wall  pressure  for  the  present  study  and 

experfm^ntard^'”'"''*  the 

In  Fig.  12  mean  s^amwise  velocity  profiles  are  compared  with  the  experimental  dau  at  several 
streamwtse  stations.  The  plots  serve  to  illustrate  the  changes  in  boundary  laver  velocity  profiles  which 
T  of^teraction  flowfield.  Fig.  12  at  x/5„=-1.6  is 'the  incom^  equilibrium 
turbulent  boundary  layer.  The  location  of  the  shock  wave  appears  as  a  kink  in  some  velocity  profiles 
downstreain  of  the  compression  comer.  The  downstream  profiles  are  seen  to  recover  rapidly  from  the 
retarding  effects  of  the  imposed  adverse  pressure  gradients.  ^  ' 


6.3.  Three-dimensional  shock  wave  turbulent  boundary  layer  interactions 

The  next  example  is  the  study  of  flowfields  of  a  three-dimensional  sharp-leadine-edeed  fin  for  swem 

fin°fa  -  M  shows  the  phTsical  domain  for  a  3-D  sha^ 

fin  (a  -  20  )  with  a  general  flowfield  structure  (Fig.  13(b))  (34).  The  inlet  boundary  conditions  and  th^ 


(•)3-D20Tin 


(b)  20*  fin  iiuenciion  flowfield  stnieoires 


vw./  v.wiiipuiwonai  aomain 


H  '-«*«  Af.=2.93.  P  =20  57kPa  r  -92  39K 

adilbauc  wall  O" *‘'P 
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corresponding  flowfield  structure  are  the  same  as  in  [351.  Here  the  free  stream  Maeh  i. 
temperature  are  =  7  onw  y  _ oa  an  ir  ».  *  sm  n4ach  number  and 

of  680kPa  and  ifl  K.' respectively:  w'i’th  the'Sr'numbe?^^^^^^^ 

thickness  8^  at  the  apex  of  the  fin  is  1.4  cm.  yielding  a  Reynolds  number  Re,  =9.8  x  lOM^^deTto 
match  the  boundary  conditions  as  used  for  the  experiments  [351  the  flowfield  behind  the  fi  '  ?  f  ^ 

as  a  flat  plate  boundary  layer  such  that  the  computed  bound  ^ 

^SfiTon  r*"  ?"  w^bounclL^'S^^^  Tc 

Adlnr^'  P  ^  downstream  exit  boundaries,  the  flow  variables  are  set  free 

Adaptively  spaced  gnd  points  are  33.  41  and  31  in  the  streamwise.  spanwise  and  vemcal  dirLw' 
J^ccuvdy.  Spearal  demaats  of  Ugendre  polynomial  degree  2  in  spa«  2  are  appUed  m The  Zto 

Fig.  14  shows  the  background  flowfield  based  on  the  eeometrie  ronfi«nire.»;«.,e  j  u 
^ndnlons  described  in  Fig.  13.  as  observed  from  rhe  from  (,1°  a^^y  “  S  a”  s”" 
of  the  hidden  poraon  are  shown.  It  is  noticed  that  the  trend  is  in  reasonable  agreement  with'the  resdd 


(a)  Density  (Kg/cn  )  . 
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(c)  Tenperature  (K)  . 
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(b)  Pressure  (Pa). 


fd)  Mach  number. 


Fig.  14.  Background  flowfield  as  observed  from 


the  front  (x  -  z  plane  and  y  -  z  plane). 
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temperature  nsT being  dis^ribuletf dong  ^he"^ fllfpraTe ‘"oTm  h°"^ 

the  shock  waves  toward  the  flat-piate  bLndary.  ^  ^  ^  number  sharply  decreasing  through 

^hown  .„  Hg.  15.  ^Thc  co'iLpnn^nTveToS/v^e  - 

vonex  stretching  occurs  toward  do^stream  with  L  ^  nght-hand  side.  Clearly,  the 

vortex  centers  close  to  the  wall.  These  physical  nhenom^aT  slip  lines  and 

stream  in  agreement  with  the  skematics  shown  in  Fig  l^h)  sigmficant  toward  down- 

rig.  16  shows  the  contours  of  vonicitv  cnmnon^n,  . 

oniaty  component  m  the  spanwise  venical  planes  (x  cos  «  -  <  planes) 


f' 


‘  *••*•***  I  Hlllllllltiiiii 

. . .  ^ 

III] . . 

'******•••  . . . 

:  :  :  :  i  ;  ;  -  : 
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Fig.  15.  (Contd) 


in  the  y  cos  a-direction.  with  each  plane  identificH  oc  a  h  ^  -n. 

toward  downstream  and  moving  upward  away  from  the  sho^k  1^.  stretching  occurs  again 

within  the  boundary  layer  close  to  the  wall.  ^  concentrated 

b-  S  c  0  contours  are  presented  at  various  locations  (a:  2Sn, 

wall  with  IB  intanal,,  increasing  toward  dol„s,„am“  exacted 

in  (241  detaiIed°sn^rafmnvlTO'”°*'h*l^^*'^'°^  ™  boundary  layer  as  discussed  in  Section  4  and 

;;“';ut:?irte"::^ereTer“  -  "»•  ^ 
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concentrated  within  the  boundary  layer  close  to  the  wail.  ^  «ipwar  away  from  the  shock.  The  growth  of  vortidty  is 


7.  Conclusions 


Based  on  the  preliminary  results  obtained  for  the  dirert  mir*...  i  i  • 
generalized  Galerkin  Legendre  polynomial  spectral  element  LthnTT  “=^‘"8  ‘“e  MEl 

for  the  initial  attempt  has  been  successfullv^achieved  Flah^  appears  that  our  onginal  goal 

been  developed  are  the  maior  factor  for  thfc,.  h-  structure  schemes  which  have 

three-dimensional  computations  of  sharp  fill  wave^IlXTen^bound‘*‘^‘^°"’  preliminary  results  for 
2-D  benchmark  problems  are  satisfact^.  ^  boundary  layer  interactions  as  well  as  the 

fluauations.  unstead^eM.'^and  tuTbule^^^  P’®’'  (1)  verification  of  3-D 

Prandtl  number,  and  turbulent  Reynolds  number  turbulent  Mach  number,  turbulent 

relammanzation.  (3)  enerl  spfamm  characte^tion  of  compressibility  effects  and 

V  ;  energy  spearum  data  versus  frequency  domain  and  complete  3-D  turbulent 
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°P"™'  “">™'  01  h-p  interaction. 
Others,  they  constitute  challenging  future  tasks  parameters,  among 

direct  numerical  simulation  for  turbulent  mm  ^  ki  0°  it  is  concluded  that  the 

element  method  appeal  be  proSg.  with  the  Legendre  polynomial  spearal 
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Appendix  A 

A.l.  Analogy  between  MEI-GGM  and  FDM-T\'D 

For  simplicity  let  us  write  (12)  in  terms  of  one-dimensional  linear  functions  and  one-dimensionai 
Jacobians  a.  i  and  c  with  all  Neumann  boundary  terms  neglected.  Integrating  (12)  by  parts  we  obtain 

^44£;r'  -ii,;:;) 

Ax  ,  , 

+  +  +  +g^,  -2g;  +  g;_,) 


Neglecting  all  diffusion  terms  and  adjusting  nodal  points  arbitrarily,  we  have 

\  r  f  rt  ^  I  1 

-  ag;:,  )  - 2At/r:;  -  ag;’') 

-r.t) +r-2) 

The  FDM-TV'D  for  the  1-D  Euler  equation  is  written  as 

dt/ /  a  r  1  1  -1 

-IT  k-  - +i  -  «/-)  -y  -  U,.,)] 


(A.l) 


(A.3) 


a  -max(0,  fl)--j(a  +  |a|)  a  =  max(0,  a)  =-j  (a  -  |a|) 

Introducing  an  implicitness  parameter  s  for  the  time  derivative  on  RHS  of  (A  3)  in  the  form 

u,^v:.sMjr 

Substituting  (A. 4)  into  (A.3)  and  assuming  that 
fl-=0  a^=a 
we  obtain 

-^=3^(A£/r*'  -AG;:;)-i||^(AGr'  -zag;:;  4.  At/;-) 

■^:^(Fr  -F-.,)  -|~(f:  -  zf;., +f:.,)  (a.5) 

nf°rr- -‘r®  that  wUh  s,  =  -s/2,  s,  =  s  Ax  iF/(a  Ax),  and  - 1  for  the  coefficient 

f  ^  ‘,u  term,  we  note  that  the  MEI-GGM  formulation  and  FDM-TVD  scheme  are  analogous;  in 

S?;  rM  "  I  u'  ?  made  above.  The  implicitness  parameters  s,  and  s,  in  the 

^Ivnn^ ?h  implicitness  paramelers  s^  and  s„ 

•  ®  concept  o  scheme,  together  with  s,  and  s,,  are  expected  to  govern  complex  physical 

^enomena  such  as  turbulent  boundary  layer  interactions  with  shock  waves,  finite  rate  chemistry, 
widely  disparate  le^h  and  time  scales,  compressibility  effects  in  high  Mach  number  flows,  etc.  Note 
that  most  of  other  FDM  schemes  (such  as  McCormack  or  Beam-Warming  schemes)  arise  with  proper 
choices  of  implicitness  parameters  in  (A.l). 
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B.l.  Integrals  of  spectral  element  interpolation  functions 

aL  =  I  <pX  d/2  a:„^  =  <P^cpl  d/2  d/2  <i>\ 

ai\=Ij:4>U/2  AZ,=lKKd/2  Az,,=l<p:Aggdn 

AZ„  =  l^<P„,^„d/2  AiZp  = 

Kkun=l<i>„,,<p„dn 

"I  +53Cr;„'^a>„.,]  + -^  ^ 


with 


<>.  =  +  i’.r.a,,.  e„„  =  a,„6,„  *  b,,,b„ 


B" 

anpn 


"L  +i3^j0i..0:, +*3c,/„0:>:^  j  +^4e„„)o:,.,}  d/i 

flU.  =/„  ^53i,,J0l.,03  +^3C„„0:..,4),J  j 

CZ.„  {-M(i,a„.  +s,b,J<pZ,<p‘  +s,c„^<pZ,<p‘  j+-^(s,d,,,,  +s,e„J<pZ,<p*  ,j  d/J 
“I  {-M(^,fl„.  +  J  +s,e,j<pZ,cp:^Jdn 

oL^„  =  l  ^^3^„)0!.*>,„ +53c„„0l*>„^,j  +5.e„.)0!.*>„„.,}  dn 

=1  +^3^„)0-.*...0a  +^3S„0„*...03.J  +^(sA,r.+s,e„J<P^,^_,<p^  .j  d/J 

^*“'*  "L  +^3^,.)0„*.>!  +^3C,,„<A-*.>!  J  +-^l5,d,,„  +5.e„,J<p^^^  dfl 

=/^  {-M(^,«,>. +53^>j0„*.>; +^3s„0„*.>:,j d^ 
=£  {-M(i.«,,.  +^3*<«)0.*.>,,,  +  J3C,;„0^*.>„„J  +^4V.)«’«*->-.,.,}  <W 

<,  =  1^  {a/  0:,.0,(f;„  +  g;j  -  ~  (a,,,  +  b,j<t>Z<p^jr^.^  +  g;,j}  d/j 
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^  I  {ax  AU--  sAK  lur'  -  C„„  Ai;:r )]  -  6,.,a„J 

+  bMl  A^::'  jx:,  dr-l<p:  [axCF:  -  G- )  -  ^(a,„  -  ,  -  G". ,)]«,.  dT 

""  in  ~  Gg„)  ~[a,r,  +  -  Gg/,)!  dfi 

+  /^{AX'^„*[-x,a.>,A£/;  '  -5j(6(„ At/;”  -c,.„At/;”)] 

+  At/:”}n,  dr -I  <PL[Ar(f-  +  G^)  -  G' j]„,  dT 

~  in  1*^  Gg(,)  -  — (a,.„  +  +  Gg„)|  dn 
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for  further  follow-up  for  the  researcher  as  well  as  practicing  engineer. 
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combustion.  This  is  due  to  Sd^v  nonreactmg  flows  or  iow-speS 

atpeaes  associated  with  finite^rawT  region  rates  for  disSent 

iayw  interactions.  C^^f^  SnSi  turbulent 

^  as  Mach,  Reynolds,  Peciet,  variables 

n<^  points  and  between  time  smm  between  adjacent 

ratable  comuntatinnsi  ..i _  P*  tnxportant  factnr*  in  _ 


n^r  •  -tteynoids,  Peciet,  and  DamlroMJ:  u  variables 

n<^  points  and  between  time  smm  between  adjacent 

ra^le  computational  schemes  2  thS  selecting 

“pw.  It  IS  toward  this  idea  that  the  nn21  indicative  of  physics  of  the 
Dependent  Muted  Explidt-LpUdt 

tJ"  paper.  Here,  eJery  n^^  ndt  ^  is  developed  in 

^que  xiumericaj  scheme  accordinff^m  •  ^®nxeiit  is  provided  with  a 
c^teriz^  by  a  total  of  six  implidtn^^^^*  situations,  as 

fi^he  ^ges  of  Mach,  calculated 

between  adjacent  nodal  points  and  hot™  ^®^’  Damicohler  numbers 
depenuent  implidtness  paramete^^n^  These  flowfield- 

conipressibie,  incompressible  visco«2^^*T  i®*^®^  flowfieid  is 

ai  the  status  of  chemic5reMt’ioMT^U“'^®“‘^’i““ar,  turbulent 

of  all  variabl®s  between  *adi^J/!i°Ti'“®’ or  fluc- 

time 

tje  fluctuation  components  “fo^l'iations 

f^«*atheNavi®r-Stokessoiutio2\n^»,^  °‘»*ained  as  the  difference 

VadeS?  standard  time  avirag  “or  the  components  calcu- 

mi^K  solutions.  If  all  imSfdtSJ  na«  *  transform  of  the 

numbers  instead  of  being  calcnlafif^?“^  parameters  are  fixed  to  certain 
practically  aU  currently  available  num^riMi*^*;,^^^®^*^  information,  then 
proposed  approach  is  believed  schemes  arise  as  spedal  cases 

simulation  of  the  supersodrSmfa«.i®  ““J  applicable  to  direct  nmSS 
higi  and  ™ ‘“>5«»tnm  SraS  IS 

D“  •<>  “mpuw  SSS,“k^“^  layer  intnraciions 
fnJ^f  P^P®^  “®  °“ly  for  coaiS  mSh  Infc’-^  T*'"®"’  numerical  examples 
turbulent  microscales.  Furthw  stSS? ™^^®  without  resoivine 
*"‘■“"0“  “  combnnrtn  «  “ 


combustion  geomet- 

shock:  2  -  bounS?  ll;  J.  7 

disk.  7  —  expansion  waves;  8  — 
layer  9  —  reattachment  shock. )  (J)  Ramiet 

1  -  ob^i“SS 

waire.  j  nonxiai  shock  system:  J _ fu«i 

uueeurs;  4  ~  fud;  5  -  Lne  holders-  fi 

—  engine  cowi;  7  —  vehicle  boundaSy.  (c) 

^  Clique 

wave,  2  --  fud  injectors:  2  —  fuel-  4 

—  eogiae  cowl;  5  —  vehicie  boundary  ’ 


xVeariy  half  a  century  has  elapsed  since 

the  digital  computer  revolutionized  com¬ 
putational  techniques  in  engineering  and 
m^hematicai  physics.  During  this  time 
finite  difference  methods  (FDM)  have 
dominated  the  field  of  computational 
flmd  dynamics  (CFD)  (1-4],  whereas  the 
opposite  is  true  for  finite  element  me¬ 
thods  (FEM)  in  solid  mechanics.  In  re- 
^t  years,  however,  the  trend  toward 
finite  element  methods  in  CFD  appears 
to  be  increasingly  favorable  [5-9].  Nu- 
n^cai  methods  for  supersonic  combus¬ 
tion  mainly  with  FDM  have  been  active¬ 
ly  pursued  since  the  late  1980’s  [10-13]. 

One  of  the  most  important  questions 
in  CFD  is  how  to  deal  with  large  gra¬ 
dients  of  the  variable  (density,  veloci¬ 
ty,  pressure,  temperature,  and  source 
terms).  Rapid  changes  of  Mach  num¬ 
bers,  Reynolds  numbers,  Peclet  num¬ 
bers,  and  Damkdhler  numbers  (if  re¬ 
acting)  between  adjacent  nodal  points 
or  elements  can  be  a  crudai  factor  in 
determining  whether  the  chosen  com¬ 
putational  scheme  wiU  succeed  or  faiL 
Furthermore,  proper  treatments  for  in¬ 
compressibility  and  compressibility,  vis¬ 
cous  and  inviadd  flows,  subsonic  and 
rapersonic  flows,  laminar  and  turbulent 
flows,  nonreacting  and  reacting  flows 
are  extremely  imponant.  These  various 
flow  properties  may  be  depicted  in  ty¬ 
pical  reacting  flow  problems  for  ramjet 

Md  scramjet  combustion  as  shown  in 
Figs.  la,6,c. 

there  be  a  general  approach  to 
satisfy  all  the  requirements  mentioned 
above?  Can  a  single  mathematical  for- 
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motion  ie^  to  most  of  the  currentiy  available  computational  schemes  both  in 
FDM  ana  FEM  as  speaai  cases?  Most  imponantly,  wiU  such  an  approach  guarantee 
accuracy  and  efficiency?  In  this  paper,  our  goal  is  to  respond  to  these  questio^ 
positively.  To  this  end,  our  approach,  known  as  the  Howfieid-Dependent  Mixed 
Exphat-Impbat  (FDMEI)  method,  is  based  on  the  foUowing  procedure  (7-9): 

(a)  Write  the  Navier-Stokes  system  of  equations  in  a  conservation  form. 

(b)  l^and  the  conservation  variable  in  Taylor  series  up  to  and  including 
the  second-order  time  derivatives  of  the  conservation  variables. 

(c)  Introduce  in  step  (b)  six  different  flowfieid-dependent  impUcitness  parameters 
wb(i  are  calculated  from  the  changes  in  Mach,  Reynolds,  Peciet,  and  Dam- 
kohler  numbers  (if  reacting)  between  nodal  points  or  local  elements. 

increments  of  the 

T  ^  form  resembles  the  impHcit 

factored  scheme  of  Beam  and  Wanning  [1],  but  much  more  rigorous. 

(e)  Step  (d)  may  be  used  either  in  FDM  or  FEM. 

procedui,  a,  danibed  above  U  capable  of  resolvmg  complex 

rvacting  iowe 

(1)  Shock  wavee  in  compreesibl.  flows  are  dependent  on  changes  in  Mach  nn- 

«av.  discontmmties  are  charactetiied  by  these  changes  in  Mach  nnmber. 

“  d«nges  in  Heynolds  num- 
her  be^  no^  points  m  FDM  and  within  local  elements  in  FEM.  Incom- 
piesstbUity  conditions  are  characterised  by  these  changes  in  Reynolds  nnmber. 

a‘  *?’"  “•  “  dianges  in  both  Mach  nmn- 

SSL  "r  Pdnts  in  FDM  and  within  locai 

ni^r  ”  " 

“  dinnges  in  Peciet  nnm- 

te^een  no^^  points  in  FDM  and  within  local  elements  in  FEM.  The 
M^Stbl  ”  “  “  characterised  by  thews  changes  in 

(6)  Renawg  flows  are  dependent  on  changes  in  Damkohler  nnmber  between 

nodal  points  in  FDM  and  within  loc^demwits  i.  r™  %t. 

s-»Ynn»r.»ieem  w..  _f  ^  eicmeuts  lu  t  EM.  The  mass  source  vs 

amvective  transfer,  mass  sonrce  vs  diffusive  transfer,  heat  sonrce  vs  convective 
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diifciv,  heac’rr,^  « 

(6)  Direct  nuxnericai  simulation  in  wu*  l  « 

out  until  turbulence  length  mici^<2a  ^e  Jo"t«nimtt  ere  curried 

can  not  be  reliable  particuiariv  fnr  >,•  without  turbulence  models 

.be  cotnpu  ^ 

variables  as  described  in  (3)  above  To^l^  ^  treating  high  gradients  of 
Legendre  polynomial  spectral  modes  mav  be  addeT*  turbulence  calculations, 
not  the  spectral  mode  approach  is  aS^t^^  .  “ 
efficiency  remains  to  bf  seen 

the  example  problems  in  this  paper  are  n  1  computer  time, 

resoiutions.  ^  ^  ^  mtended  for  DXS  microscaie 

1  followed  b, 

io  Section  4,  end  c=nclu.Lg  remiC  st^on  i  ’  ”” 

2  mathematical  FORMULATIONS 

dt  3*,-  9x,-~®’  (1) 

variables,  di^on  flux  variables,  ald^so^^t^^  variables,  convection  flux 

’  source  terms,  respectively, 


Us= 


G,= 


P 

pVj 

pE 

pYk 

0 


P,= 


-TV/Vj  -  kTi  -  ZpCpkTDk„,Yk,i 

-pEkmYk,i 


pVi 

pviVj+pSij 
pE-Vi  +  pv,- 

pnv. 


0 

Pfi 

-■fffctD*  +  pfjVj 
V>k 


where  fj  Ykfkj  is  the  body  force  Yu  is  thp  rh  i 

pom.  enmelpy,  j,  ^  “  ‘ke  Yemeni  epedee,  it  ,he 

equations  ibr  vibtationel  end  electrLie  '"““1'  kiSusivity.  AddlUoi 

hypersonics.  energies  may  be  included  in  Eq.  (1)  i 
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Expanding  tke  conservation  variables  U  in  Taylor  series  including  tbe  first  and 
second  derivatives,  we  have 


=  IT  + 


dt  ^  ~  ~2  ^  , 

where  si  and  s,  are  the  implicitness  parameters  defined  such  that 

5Un+«t  ^un 

-  >+ Sl¬ 


at 

^Un+f, 

at2 


at 
a^ir 
ata 


Sj- 


at  ’ 

5*AU’*-^‘ 


at» 


0  <  Si  <  1 , 


0  <  S2  <  1  , 


(3) 


(4) 

(5) 


wth  convection  flux  F,-  is  a  function 

of  U  and  the  diffusion  flux  G,  is  a  function  of  both  U  and  its  gradient  U ,  Thus 
we  have:  '  ’ 

dGi  „ 

_  ^7  +  ®’  (6) 

a*u 


au 

-£Ei- 

dt  ~ 

dx,- 

^±( 

dxi\ 

.'■’irj 

where  the  conve«ion  Jacobian  a,-,  the  diffiision  Jacobian  b,-,  the  diffusion  gradient 
Jacobian  c,-^,  and  the  source  Jacobian  d  are  defined  as  ^ 


af=—  b-^ 

^  au  ’  ’  ”  W 


Cij  = 


dGi 

au“ 


d  =  ^ 

au 


(8) 


fuming  the  product  of  the  diffusion  gradient  Jacobian  with  third  order  spatial 
denvatives  to  be  negligible,  it  follows  from  Eqs.  (1)  through  (7)  that 


aAG?'^* 

ax,-  axi 


~  AP"t- 


iV 


+ ,  +  b.)  j 

~  — ^')]|+o{ai’).  (s) 

In  to  prwdo  different  impUatnere  (different  nmnericni  trentments  or  schemes) 
to  different  physicel  quantities,  we  reessisn  s,  and  s,  associated  with  the  diffusion 
and  source  terms,  respectively, 


Si  AG,  saAG, , 

AG,  =>  a,  AG, , 


3i  AB  =>  dsAB , 


S}AB  =>  SflAB , 


(10) 


with  the  various  implicitness  parameters  defined  as 
3i  s  first  order  convection  implicitness  parameter, 

S2  =  second  order  convection  impiidtnesa  parameter, 

33  *  first  order  diffhsion  implicitness  parameter, 

34  =  second  order  diffusion  impUdtness  parameter, 

35  =  first  order  source  term  impiidtness  parameter, 

3«  *  second  order  source  term  impiidtness  parameter. 

den»d.!r‘  parameters  Sj,  S3,  and  s,  wiU  be  shown  to  be  fiowfield 

rs  :£..r  =  s 


V  *1- 


2  rl  - — 


+  ^4 


r  /  5’fa,b,-  +  b,-b^)AU"+»  \ 

J 

J  /  ab,AU"+>  d»c.vAir'+‘  \  1 

[\  dxidxj  j 

—  Q  t  ] 

\  dxidxj  J 

^ - dAir+‘ 


■1) 


with 


Ux,.  ax,.  BJ-  — + 
,fSF^  dQy>  M 


(ll) 


Raairanging  Eq.  (u),  we  obtain 

d 


=  ^Air+'  =.  dAtr*' , 


or 


R  =  AAIT-'  +  ^  {EiAU"+>)  +  (E,.jAir-‘)  +  q-  +  0(^3) 

=  -Q" . 


(12) 


(13) 
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with 


Af2 

A  —  I  +  sgd , 

Et  =  Ai(3ia,-  -h  asb,-)  +  -^;36d(a,-  +  b,)  +  ajda,-  +  54db,'] , 

T.  r 

E,-,-  =  AtSidj  -  —  la2(a,a,-  +  b,aj)  +  34(a,b;  +  b,b^  -  ddj)]  , 


(14) 

(15) 

(16) 


^  dZi 


•d  )  (17  +  G?)  +  ^(a,-  +  bOB" 


(^1+^ 

5^  +  'xXU  -  s?)]  -  (a<  +  B- . 


(17) 


An  alternative  scheme  is  to  allow  the  source  term  in  the  left  hand  side  (LHS) 
Eq.  (13)  to  lag  from  n  +  1  to  n  so  that  Eq.  (13)  may  be  written 


01 


as 


with 


^  +  «?)  +  ^<X  +  bOB' 


(18) 


a* 


dzidxj 

/a. 

-  I  Atas  +  ”^36 


At* 

— {«.  +  bO(F74-G7) 


(‘ 


Af2 


dAir-  lAt  +  i^djB". 


(19) 


Warming  scheme  [1]  can  be  written  in  the  form  identical  to 
hq.  (IS)  with  the  following  definitions  of  E,-,  E,-,-  and  Q": 


E,‘  —  mAt(a,-  +  b,-) ,  with  m  - - 

1  +  ^’ 

Eij  =  mditdj , 

q»  =  AL/'£!I^a^V. 


(20) 

(21) 

(22) 

are 


wh«e  the  cross  derivative  terms  appearing  in  Q"  for  the  Beam- Wanning  scheme 
indudeQ  in  the  second  derivative  terms  on  the  LHS.  The  Beam- Wanning  scheme  is 
seen  to  be  a  special  case  of  the  FDMEI  equations  if  we  set  si  =  53  =  m,  =  34  = 
35  _  34  =  0,  in  Eq.  (18),  with  adjustments  of  Q"  on  the  right  hand  side  (RHS)  as 


AN..  ComputaMonai  Appto^j,  .i.t  Fl,.<irid.Dep.n„.n.  ImplicitM^. 

u  c.  .  ...  ..  FBMB:“SS“.=rEr^3,TL^ 

•  ,  capable  o2_  produdng  practicaUy  aU 

aT  f^>^=t/ioo'^  ' existing  FDM  and  FEM  schemes.  Some 
^  examples  are  shown  in  [9]. 

/m  Beam-Waiming 

i  /^/^<  y^/m\  approach  is  to^ 

§  OJQ  {//y  impliatness  parameters  from 

J|  2-  //  '  •////  /  I  current  flowfield  variables  at  each 

^  I  0.25  y  3 wV  ff P°“^  than  fixing 

1 1  y y  the  impliatness  parameters  to  certain 

\y  ^yyyy  ••..  J  P*’«i«enmned  numbers  and  using  them 

0  025  OJ  ■  '0  75  I  ®^“^®5°^*^°“^^’^®Pcctiveof 

0.639\q  fi°^eid  variation  from  one  po- 

Fust  order  impUcitmas  parameters  j,.  another,  in  FDMEI  scheme, these 

implicitness  parameters  may  be  deter- 

Pimm^  J  B-._-  m«..j  r__ _ .  ,  .  ^  «eier- 


J  0.75 

s 

I OJO 

|5 

1 1 
S'?. 


‘1“  fet  temporal  base. 

-»=  “  ftg-  2.  The  Saai  wlaes 

..•maijl-,,,  l.).a,  =„awi_4’f ,,  ^  unpiiatness  paiame.eis  a.  any  po- 

=  t  “  >>'  obtained  a. 

<  1.  D  —  fleem-Wanmng  ‘'“Ay'rage  of  both  spatial  and  .empoial 

contri  fautions: 

Conveetion  Imoiicitness 


1)  r>a 

0 

^  ^min  =  0 


■»a  =  «i ,  0  <  n  <  1 


r  s  y^*^***  ^mtn 

,  ’  (24) 

wnere  tae  maximum  and  miTiiTnst'pn  MaM,  nnTTsj, 

nodal  points  in  FDM  or  within  a  iocaJ  fi  . .  calculated  between  adjacent 

parameters  and  between  the  time  sten  J  *P«^al  implicitness 

parmneten  and  a  is  a  nset-spedaed  sLu  rarber’'(a"-  n^irT”* 
ai  IS  directly  related  to  the  flowfield  who«.„  J  ^ 

Thepnmatyroieofs,  is  M,  ens^T.  “  '>  "a  = 

theconrection  giadienB,  whereas  that  of  s  '=!'  Pfopariy  accommodating 

Stability.  ^^ficiai  viscosity,  for  solution 


Diffusion  fmoiicitness  Pammpr^^,.. 


S3  = 


s  >  0 

^  <  0,  J2e„in  #  0 , 
—  0  , 


with 


or 


or  ^  0 

-^^min  —  0 


T^O  54  =  55,0<n<l 


(25) 


(28a) 


Peclet  numbers  ire  "^ters  or  mammnm  ind  minimnin 

parameters,  ajid  /?  «  a  user-specifie/.,™  V'  ^  ^®“Porai  impUcitness 

ir.  iirge,  i.  poS^pS  -  MD-  H  ■enfpem.^ 

*««e  tbe  diSnsion  impildtnes.  pi^Z  of  Reynolds  numbers 

el>oeen,  is  obtiined  eitber  from  E,.  (26^(26M  ?*.  ‘T' 
f^mog  rbe  solution  iccnnny  bv  tLfJZ  j  ^'  >*““««  e.  =  sS  wHb  s, 

^  s,  pliys  rb.  SiibSr^ir” 

U?*?  Perm  - - 

- .«  n.  c.™..  ^ 

f  itlin/' #  1  ^  4  . 


with 


nuii(t,l)  i>y 

*•■=(  0  (<T.a„n, 

^  -Oarnm  =  0 


•*«  —  "S?  1 0  <  ra  <  1 


(27 


,  „  -  Qe^., 

•kere  , be  minimum  ind  minimum  Di^r.  ,. 

“  “.^eok  es  for  spitial  end  tempoml  imitiStn™*''”  similirily 

speaned  smail  number  (7^0  ni^  P^rameters^  and  7  is  a  user- 

•koe.  ibr  convection  ^<> 

^ntrolling  the  solution  with 

average  of  both  spatial  and  tempor^^J2dtnl/n  ""  ^^biUty,  respectively.  The 
m  ‘^utation.  at  any  point  (elementT^d  Sme" 
aalationships  between  all  ohvaiMi  iu 

tremmeuK  ir.  cbinmerized  by  “^P-otog  -ounericnl 

0—  (..  es.  ss,  md  J seconSt^dS  L^n^ 


vi«osi.y  m  ordar  .o  praa  Ja  .ha  t  “““'  »f  “>«=rtaai 

for  tile  second  order  impUdiness  .  ^carzcy.  Note  that  the  definitions 

in  (7.  8]  in  order  to  me^t  1  "  ntodihed  front  those  reported 

second  order  impHcitness  parameters  should  be  th^!? 

the  first  order  implidtness  parameters  fs,  =  i  _  compliances  of 

such  that  the  second  order  implidtness  oarametPras^aV  ^  *  1  “  «$)  [8] 

respectively,  for  the  minimum  and  maximum  v»l«  ““““m  and  minimum, 
parameters.  Unfortunately,  such  definition  r  it  ^  implidtness 

for  U..  high  vaiaas  of  ,l"’ ST orte “oSS  “  ““ 

l™tiagraiaaa(OJ)ofthas«condortarimpEdtMs’’n^™“”'  <1>« 

that  sj  =  maxfl  —  Q  st  «#/.  --  •  ^  parameters  were  provided  such 

both  first  and  second  order  parame^s^^Srls'^  However,  it  was  noted  that 
««ramaa  a.  aaro  aad  uni.,  wirtlTa^^d  o^^*'  ?  '>'>“■ 

«»oaahiy  larg.  for  ail  Jaaa  of  ,h.  ^ord^JS?^r““  '“'‘“S 

^  sacond  order  impHcitaass  paramelara  dTan  S,  P“»matars.  Tima, 

tatttiona  of  the  Sm  order  impiidtnasa  coatimooa 

The  range  of  the  constant  n  is  0^<  n  <  i  satisfying  these  requirements. 

optimum,  exhibiting  the  best  converEence^^°T  ”  “  *be 

the  aiafflpla  problaL  p  teaBnad  ip  ^opT  ”  “S*"  CH.  pupibara  ip 

ollowipg  varioua  PPmeiS's™mM'^be^''''”^‘*'*^  oapable  of 

aa  fbUoara;  «  >>'  »otop.pt.cany  gaparatad,  aa  spptmarizad 

(1)  The  first  order  implicitnp^^  par3jnntm-«  «  ^ 

nomena  such  as  shock  waves  and  turbulLM  Th«r”^  ^ 

from  the  changes  of  local  Mach  nnmK  .i  n  Parameters  as  calculated 

within  each  elLent  anTare^dic^'^, 

The  contours  of  these  parameters  d  i  ^  element  flowfieids. 

with  both  Si  and  S3  bdne  ian^T?  resemble  the  fiowfieids  themselves, 

Pda  of  a,  apd  a,  is  to  previda  coptpptatiopd  aS^^  *” 

(2)  The  second  order  implldtups*  . 

ent.  However,  their  primary  role  S4  are  also  fiowfieid  depend- 

(^fidai  viscosity)  as  they  were 

time  derivative  term  of  the  Tavior  introduced  into  the  second  order 

variables  IP^i.  ^le 

stability.  ^  provide  computational 

o’it'pUS  ^ » «■“ 

au.  iha  cotppptatiopai  schapia  is  aptomaticaiiy  aitarad  to 
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take  this  effect  iato  account,  with  the  governing  equations  being  predonunantiv 
paraoolic-^iliptic.  Note  that  these  effects  are  confined  at  U’*+\  not  at  U". 

(4)  The  53  terms  «e  associated  with  diffusion.  Thus,  with  S3  S  0,  the  effect  of 

f  P“smet«s  ..  and  a,  are  nonrero,  this  indicates 

a  typical  sunatjon  for  the  misted  hyperboUc,  paraboUc  and  elliptic  natnmTf 

^  baeter-Stokes  system  of  epnations,  with  convection  and  i^ffhsion  blv 
^«lly  important.  This  is  the  case  for  incompressible  ho,™  at  W  soe^ 
^  nmqne  property  of  the  FDMEI  scheme  is  its  canability  to  control  prSre 
Rations  adepnately  without  resorting  to  the  separate  hyperbolic  effiutic 
p^uie  equation  lor  pressure  corrections.  The  capability  ofTOMn  sthm. 
»  h^dle  in^mpressib  e  Sows  is  acbiev«i  by  a  diicate  bali™ 

^bL  ffrlTr  (or  Peclet 

rE,.T4  s  “““Ptaaaiblo  (M  =  0).  the  crtU  ptu 

j  whereas  the  implicitness  parameter  st  is  tn  Ko 

between  converting  anri  rms.  These  terms  allow  interactions 

cd“lnrr  ■“  ““  -™«ble  and/or  viscous 

Such  ui  sp^' ““P- 

the  waU.  ^  temperature  compressible  flows  close  to 

which  are  fnnctteiri? 

term  implicitness  narametPr  •  ^  numbers,  the  first  order  source 

the  resulting  equations  to  obtain  “  *^“^6  with  the  stiffiiess  of 

other  hand,  the  second  order  so^TeTS'u  •“  “« 

to  the  stability  of  solution?  rt  •  impiiatuess  parameter  sg  contribute 

wiU  adjust  the  reaction  rate  terlL'ta”  ‘^V*" 

time  to  the  reaction  time  in  S  ratt^lT 

solatioas  with  computational  «^i™  'P' 
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(8)  Various  definitions  of  Pedet  number  and  Damkohier  numbers  (Table  1)  bet¬ 
ween  the  energy  and  spedes  equations  snouid  be  checked.  Whichever  definition 
provides  larger  values  of  53  and  S5  must.be  used. 

(9)  The  transition  to  turbulence  is  a  natural  flow  process  as  the  Reynolds  number 
increases,  causing  the  gradients  of  any  or  ail  flow  variables  to  increase.  This 
phenomenon  is  the  physical  instabiUty  and  is  detected  by  the  increase  of  s,  if 
tteflow  is  incompresdble,  but  by  both  53  and  sj  if  the  flow  is  compressible. 
Such  physical  instability  is  likely  to  trigger  the  numerical  instability,  but 
will  be  coxmtered  by  the  second  order  impiidtness  parameters  33  and/or  sa 
to  ensure  numerical  stability  automatically.  In  this  process,  these  flowfieid 
dependent  impUdtness  parameters  are  capable  of  capturing  rdaminaxization 
compressimlity  effect  or  diiatational  turbulent  energy  dissipation,  and  turbu- 
lent  unsteady  fluctuations. 

(10)  important  contribution  of  the  first  order  impUdtness  parameters  is  the  fact 
that  they  can  be  used  as  error  indicators  for  adaptive  mesh  generations.  That 
IS,  the  larger  the  impUdtness  parameters  the  higher  the  gradients  of  any  flow 
gabies.  Whichever  governs  (largest  first  order  impUdtness  parameters)  wiU 
locate  the  need  for  mesh  refinements.  In  this  case,  aU  variables  (density 
w^aty,  pressure,  temperature,  spedes  mass  fraction)  partidpate  in  resolving 
the  adaptive  mesh,  contrary  to  the  conventional  definitions  of  the  error 


(11)  Physi^y,  the  impUdtness  parameters  wiU  influence  the  magnitudes  of  Jacobi¬ 
an*.  Thus,  Item  8  above  may  be  modified  so  that  the  diffusion  impUdtness  pa¬ 
rameters  33  and  34  as  calculated  from  Reynolds  number  and  Pedet  number  can 
be  appaed  to  the  Jacobians  (a,-,  b<,  c,-,).  corresponding  to  the  momentum  equa- 
lons  aad  eneip  equation,  respectively.  Furthermore,  two  different  definitions 
of  Pedet  number  (Per,  Pt„)  would  require  the  33  and  34  as  calculated  from 
the  ener^  and  spedes  equations  to  be  appUed  to  the  corresponding  terms  of 

T'  for  the  source  term  impUdtness  parameters 

3,  and  3e  should  be  foUowea  for  the  source  term  Jacobian  d.  In  this  way,  high 
temperature  gr^ents  arising  from  the  momentum  and  energy  equations  and 

the  tote  rate  chemistry  governed  by  the  energy  and  spedes  equations  can  be 
resolved  accordingly. 

^  “"-y  FDM  or 

.  Galerkin  approximations  of  FEM  lead  to  the  results  of 

central  differences  of  FDM.  However,  the  main  difference  between  FDM  and  FEM 
Mses  wnen  integration  by  parts  is  performed  in  FEM  and  the  expUdt  terms  of 
Neimann  noundary  conditions  “naturaUy”  appear  as  boundary  integral  forms.  Thus, 
eumaan  boundary  conditions  can  be  directly  spedfied  at  boundaries  in  FEM. 
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_ TMe  1.  Definitions  of  nondimentionai  flow  variables 


p{v  •  V)v  =  - Vp  +  n[V^v  +  r VfV  •  v)! 

-  o 

A  B  - - - - - 

C 

^  /  T  \ 

v.p,|c^ir-v^-v.  t)Dvir/[c^dr\=.f'Htv,i, 

To  F  \  I 


V .  (pYkv)  -  V .  {pDVYk)  =  Wk 


Mach  number 

Af 

Be3mQids  number 

Re 

Peciet  number,  I 

Per 

Peciet  number,  II 

Pen 

Damicohler  number,  I 

Eat 

Damkohler  number,  II 

Dart 

Damkohler  number,  HI 

Earn 

Damkohler  number,  IV 

Dajv 

Damkohler  number,  V 

Day 

M  C  viKotis  force 


E  convective  heat  transfer 


F  conductive  heat  transfer 


convective-  mass  transfer 


diffusive  mass  transfer 


1  ^ _ mass  source 

mtYk 


pDYk 


co&vective  transfer 


mass  source 


difiusive  transfer 


convective  heat  transfer 


^  conductive  heat  transfer 


^  _  heat  source 


qL^  N 


HD 


wug  nowneld-Dependent  Impiiciineas  Algorithm 
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cumbersome  process  must  be  taken 

bo^dary  layer  imliSti^sw^eMmp^^^  ^  shock  wave  turbulent 

and  invisdd  hows,  and  laminar  and  fl  ^compressible  flows,  viscous 

the  fiowfield  domain,  a  computational  scheS' inteLd2^*^^*T^®^  throughout 
physics  and  that  does  not  account  for  ot^^r  ‘  0“  type  of  flow 

For  example,  the  flow  dose  to  the  waU  in  sho?“  Phenomena  wiU  faU. 

interaaions  is  incompressible  (M  <  0  1)  whereaT^'^  turbulent  boundary  layer 
is  compressible  (supersonic  or  ^pej^’ 

mvisaa  flows.  In  between  these  two  extremos  th  ‘diange  to 

osc^ating  back  and  forth  across  the  boundary  lave^^^f  continuously, 

leading  edge  and  bow  shocks.  At  any  eiven^J«^  ^  ^  entropy,  and 

gradients  of  each  variable  (density,  pressure  ®ienient. 

small  or  very  large,  so  large  that  practic^v  t  i  .  and  temperature)  may  be  very 
methods  may  fail,  fo  order  to  succeed  it  is  npr  available  computational 

everywhere  be  identified  and  so  recoded  with^eJJf  Physics 

accorded  to  each  and  every  comnutatif^^  “’7^  ^  ^  computational  schemes 
such  accommodations  are  available  in  F^MEI^^^e! 

3  implementation  and  computational  process 
Stokes  system  of ’equati^^™2fiS.ne°Monyt*  Taylor  series-modified  Navier- 

finite  volume  method  (FVMl  For  FDM  or  to  the 

spatial  derivatives  order  and  second  order 

may  write  for  any  variable  u  as  ^  difference  schemes.  For  example,  we 

S  l«=  +  »■■-. 

with  analogous  formulas  for  the  « 

ext^ons  to  three-dimensional  problems  two-dimensions  or  corresponding 

For  applications  to  FEM  we  beein  *K 

variables  and  source  terms  as  a  Une^roJk^^^^'^^  conservation  and  flux 
nodal  values  of  these  variables.  mbination  of  trial  functions  with  the 

Applying  the  generalized  Galerkin  annr«  •  •  ' 

oalerkm  approximations  to  Eq.  (13)  we  obtain: 


(29) 


(31) 


X*“^U,r„G,-,B)dn  =  o 


(32) 
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or 


where 


~  la  'Iri'SS  (1  +  AtSs  +  ^-S6)^r»  , 

B^dr,  =  -|A«(3ia,v,  +  Sabirs)  T  ^[s2drtCiu  +  36<irt(a,w  T  6,t,) 

.  +  a4<ir«6tti]|^a.«^iJ  -  |Aii3C,-j„  ■  — [a2(a,>fajt,  -f  6,vtaj„) 

+  34(a,„ijt4  +  dCl 

+  «36,>.)  T  •~i32<irta,u  +  SedrtiOiu  T  bit,) 
+  a4<irt6«.]|  iaifl  +|Ats3C,-,>,  -  ^{a2(o,v,aj„  + 

+  uioirtbjtt  +  ifrtijt,  -  tirtCiju)]!  iaia. 


(33) 

(34) 


(3S) 


,  1 

^ 

2  (•^i» 


A<2 


♦fll  dft. 


(36) 


^  =  /r  {  h‘  +  ^,.) 


A^^ 


§a$/J 


!  '  t  .  -  1 

■'■  2  ■*■  *’■'■•)  (^i»  +  <53j.)  iaiaj  I  nj  <ir , 


(37) 


where  aJi  Jacobiajis  must  be  UDdatPfi  at  *4.  ■ 

Neumann  boundary  trial  and  test  funeti  •  ^  represents  the 

number  and  r,  the  Siobai  node 

For  three  dim^oL  7, fi  12  3  S^oL 

identification  of  each  Jacofaia^  matrix  (a,  a'”?  imply  directional 

n22»  C23,  cai,  C32,  C33)  with  r,  s  =  1  2  3  4  ‘^i3»  ^21, 

Jacobian  matrices.  Note  also  that  f  the  5  x  5 

source  terms  to  the  RHS  of  Eq.  (33)  instead  of  1  move  the  delta 

step  behind.  These  indices  can  be  reduced  •  -i  ^^®  ®°°^ce  terms  one  time 

given  by  Eq.  (33)  are  referred  to  as  FDMEI-rSl^h  e^l^iatious 

-I't  ‘k' 

integrals.  It  is  parti,Saxiy  ad^t  ‘L3  fha^  N  ^ 
through  rensvaluation  of  Jacobians  n^al  to  the^”^  conditions 

be  added  to  the  boundary  nodes  for  the  boundary  surfaces  can  simply 

th.  o.iar  hand,  ail  N«SL 

sonrca  tema.  Thew  featniea  .i,  ai^  in  “ 

boundary  conditions  can  be  handled  bv  devi^f^’  ^  implementations  of  Neumann 
boundary  nodes.  ^  vising  speaal  forms  of  finite  differences  at 

Similar  results  are  obtained  either  by  FDM  or  FP“U  «»«>(. 
tions  aerived  primarily  from  the  FDM^  Eos  r  ^  ^  accu^  of  computa- 

Rnynmds  number  (say  around  Re  >  1Q»^  ^°^ever,  with  the  increase  of 

with  appUcations  of  special  Wio^  3^!“  --X  -crease 

modes  cnaracterizing  extremely  smaU  tiirK«i  polynomials  of  high  degree 

high  frequency  modes  can  be  achieved  bv  ^^^*““^-*1011  ol  such 

nodes  of  isoparametric  finite  elmms  AdZnT^?  “nier 

as  n^ed  for  the  resolution  of  turbui^f 

implicitness  parameter  S3  will  play  a  cmciaJ  f  diffusion 

degrees  of  Legendre  polynomial  Ae  detennining  the  required 

superimposed  onto  isoparametric  elements  h  ^®5®““®  P°^y“o-ial  spectral  modes 

■nmbm  (iaco.^ ^1  !f  ^ 

pressure  oscillations.  This  adjustment  is  adjust  itself  to  prevent 

•mployed  for  incompressible  flows.  Otherwis®7'7hr  EDME^’^h 
shock  wave  resolutions  at  hieh  Mach  «  k  ’  *  '  DMEI  scheme  is  capable  of 

d-Uns  with  interaction,  he, particnlnrly  wdl  .nittd  for 
regions  of  high  and  low  Mach  ^  bonndary  layers  where 

case,  the  invisdd  and  viscous  intetlacr-  ®^J-olds  numbers  coexist.  In  this 
end  the  second  order  imph“me7s  !^‘‘°“^  "^®  ‘"1“  Pl-*.  To  this 

P  eters  play  the  role  of  artifidal  viscosity 
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neeaed  for  shock  wave  resolutions  in  the  presence  of  how  diffnsion  due  to  physical 
viscosity. 

In  order  to  understand  how  the  FDMEI  scheme  handies  computations  involving 
both  compressible  and  incompressible  flows  fundamental  definitions  of  pressure  must 
be  recognized.  Consider  in  the  following  that  the  fitiid  is  a  perfect  gas  and  that  the 
total  energy  is  given  by 

E-Cj,T-^  +  ^ViVi.  (38) 

The  momentum  equation  for  steady  state  incompressible  rotational  flow  may  be 
integrated  to  give 

-  j[t^  dxi , 

1 

p  +  -pvjVj  =po  +  W,  (39) 

with 

“  m  dzi , 

where  uk  is  the  component  of  a  vorticity  vector,  po  is  the  constant  of  integration, 
and  m  denotes  the  spatial  dimension.  Combining  Eqs.  (40)  and  (41)  leads  to  the 
following  relationship: 

Po  =  p(cpr  +  V,-Vi  -  E)-W.  (40) 

If  po  ^  given  fay  Eq.  (40)  remains  a  constant,  equivalent  to  a  stagnation  pressure, 
then  the  compressible  flow  as  assumed  in  the  conservation  form  of  the  Navier- 
Stokes  system  of  equations  has  now  been  turned  into  an  incompressible  flow,  which 
is  expected  to  occur  when  the  flow  velocity  is  sufficiently  reduced  (approximately 
0«1  <  Af  <  0.3  for  air).  Thus,  Eq.  (40)  may  serve  as  an  equivalent  equation  of  state 
for  an  incompressible  flow.  This  can  be  identifled  element  by  element  for  the  entire 
domain.  Note  that  conservation  of  mass  is  achieved  for  incompressible  flows  with  p© 
in  Eq.  (40)  being  constant,  thus  keeping  the  pressure  from  oscillating. 

Once  the  Navier-Stokes  solution  via  FDMEI  Is  carried  out  and  all  flow  variables 
determined,  then  we  compute  fluctuations  of  any  variable  /, 

/  =  (41) 

where  /  and  Y  denote  the  Navier-Stokes  solution  and  its  time  average,  respectively. 
This  process  may  be  replaced  by  the  fast  Fourier  transform  of  the  Navier-Stokes 
solution.  Unsteady  turbulence  statistics  (turbulent  kinetic  energy,  Reynolds  stresses, 
and  various  energy  spectra)  can  be  calculated  once  the  fluctuation  quantities  of  all 
variables  are  determined. 

Before  we  demonstrate  numerical  examples,  let  us  sumznaiize  why  the  FDMFI 
scheme  is  capable  of  handling  low  speed  and  high  speed  and*  compressible  and 
incompressible  flows,  including  shock  waves  and  turbulent  flows: 
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(1)  How  is  the  transition  ironi  incompressible  how  to  compressible  how  naturailv 
and  automatically  accommodated  without  using  two  separate  equations  or 
two  separate  codes?  This  process  is  dictated  by  the  first  order  convection 
impUcitness  parameter  si  as  refiected  by  the  Mach  number  changes  and  the 
expression  of  the  stagnation  pressure. 

(2)  How  is  the  shock  wave  captured?  As  the  Mach  number  increases  and  its 
discontinuity  is  abrupt,  the  s,  terms  associated  with  second  order  derivatives 
together  with  squares  of  the  convection  Jacobian  provide  adequate  numerical 
viscosities  through  second  order  derivatives,  similarly  as  the  Lax-Wendroff 
scheme. 


(3)  How  is  the  transition  from  iaminar  to  turbulent  hows  naturally  and  auto- 

process  is  governed  by  the  first  and  second 
orter  diffusion  impUcitness  parameters  (sj  and  s^)  as  calculated  from  the 
changes  of  the  Reynolds  number.  The  terms  associated  with  S3  and  S4  are 
responsible  for  fluctuations  of  velocities,  with  the  values  of  these  impUcitness 
Pieters  increasing  with  intensities  of  turbulence  in  conjuction  with  the 
diffusion  gradient  Jacobian  and  the  squares  of  the  diffusion  Jacobian.  This 
process  aUows  the  Navier-Stokes  solutions  to  contain  fluctuations  which  can 
be  extracted  by  subtracting  the  time  averages  of  the  Navier-Stokes  solutions. 


(4)  How  do  the  interactions  between  convection  and  diffusion  take  place?  Changes 
of  numbers  and  Reynolds  numbers  as  reflected  by  both  convection 

and  diffusion  impUcitness  parameters  dose  to  the  waU  contribute  to  the 
unsteadiness.  Away  from  the  waU,  they  contribute  to  the  transition  between 
incompressible  to  compressible  flows. 


4  APPLICATIONS 

We  examine  first  a  nonreacting  flow  problem  with  shock  wave  turbulent  boundary 

f°Howed  by  the  reacting  flow  with  a  flat 


4.1  Saperaonic  Nonreacting  Flow  on  a  Compression  Corner 

demonstrate  cdcuiations  of  supersonic  flow  on  a  compression 
O  boundary  conditions  (nondimensionalixed)  are  p  =  1,  M  =  2.25, 
~  ^  adiabatic  waU  condition.  The 

tTI  mean  flowfieids.for  the  comopression  comer  are  shown  in 

ig.  3a.  In  these  calculations,  ail  perturbation(fluctuation)  variables  are  determined 
from  time  averages  of  the  Navier-Stokes  solutions  according  to  Eq.  (41). 
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02m  “  - - 


.10256  m 


Prasure 


Tmmperartn 


Machmunber 


geometry  and  flowfielda°°(Tf^uctLt^^  comer,  (a)  Compression  comer 

r  *0.102  m,  y  =  0.004  m  ^  fluctuation  velocities,  i  -  *  =  0.102  m,  y  =  0.001  m:  2  - 


The  horizontal  and  venicai  perturbation  velocities  tu'  and  f/)  • 

to  the  wail  (x  =  0  10256  m  u  —  n  nm  m  j  ^  locations  dose 

„  -  n  na  ^  V  -  0.001  mj  and  away  from  the  waU  (i  =  n  iQo^f?  m 

y  —  0.04  m)  arG  shown  m  Fi^  vu.**  /  •  .  w-iu^oo  m, 

«gmau..ly  l«s  dole  «  .he  weil. 

s.e^^^.^  Tleee  ..ends  e.e  .e«ec.ed  in  .he  .nshnlenc. 


(I  ““®'  '‘“"il’n.ions  i.  .he  locations  opstream  of  the  come. 

(X  -  0.0513  m)  and  downstream  of  the  comer  fr  -  n  ittf  \  ^  comer 

observe  that  the  turbulent  kinetic  enerirv  «  ’  shown  in  Fig.  5.  We 

larger  than  the  upstream  Xo  turhul  ^1^®  comer  is  considerabiv 

frequences  us.  power  soectral  denshv^  calculations  (wave  numbers  or 

microscaies  are  not  resoivied  in  this  exaLpk.^"^"'^'®'^  “  turbulence 


Figure  i.  Supersonic  nonreactinv  flow  on 
compression  corner.  Reynoida  sir^? 


It  should  be  noted  that  the  above  re¬ 
sults  obtained  without  turbulence  mo- 
dels  or  without  the  standard  DNS  solu¬ 
tions  (neither  spectral  nor  DNS  mesh  re¬ 
finements)  are  regarded  as  the  consequ- 
Mce  of  the  tim6-averaging.of  the  FDMEI 
Navier-Stokes  solutions.  This  imphes 
that  the  fluctuation  of  variables  between 
nodal  points  and  between  time  steps  as 
reflected  in  terms  of  the  imphdtness  pa- 
rmeters  (s.)  have  contributed  to  these 
physical  phenomena,  with  compressibi¬ 
lity  and  shock  waves  dictated  by  the 
Mach  number-dependent  si,  and  with 
incompressibility  and  turbulent  fluctua¬ 
tion  dictated  by  the  Reynolds  number 
or  Pedet  number-dependent  S3.  The  si¬ 
multaneous  partidpation  of  and  S3 
are  also  responsible  for  shock  wave  tur¬ 
bulent  bound^  layer  interactions. 

results  of  the 
FDMEI  scheme  with  the  k-E  turbulent 
modd  and  experimantal  data  is  shown 
m  Fig.  6.  It  is  seen  that  the  FDMEI 
results  compare  more  favorably  with 
those  of  measurements  [14]. 
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Figure  5.  Supersonic  nonreacting  flow  on  a 
compression  corner.  Turbulent  kmetic  ener¬ 
gy- i  —  r  =  0.051  m,  2— rs=  0.1333  m 


4.2  Supersonic  Reacting  Flow 
with  Transverse  Fuel  Injector 

We  consider  an  example  of  transverse 
hydrogen  fuel  injector  with  9  species  and 
18  reaction  equations  and  a  rectanglar 
geometry  of  3615  isoparametric  finite 
elements  studied  in  [8].  The  primary 
air  flow  is  set  at  =  1,  Poo  = 
0.1  MPa,  and  =  1000  K.  The  s^on- 
dary  hydrogen  jet  transversely  injected 
through  a  slot  of  2  mm  is  provided  at 
M  =  1.  P  =:  0.2  MPa,  and  T  =  300  K. 
Initially,  the  frozen  flow  with  the  mass 
fractions  of  0.095  Oj  and  0.905  N2  is 
analyzed,  which  is  then  used  for  reacting 
flows. 


Figures  7a, fi  show  density  contours  for  the  frozen  and  reacting  cases.  Notice  that 
the  density  ch^ges  steadily  toward  vertical  directions  downstream  of  the  slot  for  the 
^  of  froz^  flow.  However,  for  the  reacting  flow,  there  is  a  considerable  undulation 
due  to  the  formation  of  production  species  and  the  effect  of  dilation  and  dissipation 


In  Fgs.  7c, d  the  pressure  and  temperature  variations  for  the  reacting  flow  with 
an  isothermal  w^  (300  K)  are  shown.  It  is  evident  that  the  upstream  boundary  layer 
separation  ^d  bow  shocks,  Mach  disk  and  barrel  shock,  downstream  boun^  layer’ 
and  reattachment  shocks  appear  to  be  smeared  due  to  chemical  reactions.  Expansion’ 
v^ves  appew  downstream  of  the  injection  slot  as  the  boundary  layer  is  reattached 
(  ig.  7c).  The  stattc  temperature  increases  to  the  maximum  below  the  upstream 
bounuary  layer  slightly  toward  upstream  of  the  injection  slot  as  a  result  of  of  the 
^he^c  re^tions  (Fig.  7d).  Downstream  of  the  injector,  cold  unbumed  gases 
flowng  out  of  the  injector  reduce  the  static  temperature  drastically.  There  is  an 
evidence  of  sudden  increase  of  temperature  downstream  of  the  reattachment  shock 

*°Sether  with  an  increased  oxygen  concentration 


n  the  main  products 

OH  anu  HjO  are  shown  in  Figs.  7e-h.  The  large  scale  recirculation  upstream  of 

the  injection  transports  the  injected  hydrogen  within  the  upstream  boundary  layer 
J^ation  region,  resulting  in  a  significant  amount  of  OH  and  HjO  species  therein, 
^ese  main  species  togeth«  with  intermediate  species  are  carried  downstream  of 

the  rearcuiation  region.  The  exothermic  energy  thus  created  within  the  miving 
layer  contnfautes  to  the  thrust  force.  uiijuub 


A  New  Computationai  Approach  with  Flowfieid- Dependent  Impiicitpeas  Algorithm  487 


0.081-  ^  2 
0  i 


g  0.06  - 

I 

%  0.041- 

u 

0.02- 


0.08  -  ‘ 

i 

0.06 


.§  I 

■2  0.041- 

lacL 


•02  0  0.2  0.4  0.6  0.8  I 
UAJ^ - 


-02  0  02  0.4  0.6  0.8 

UAJ„ - 


5  CONCLUDING  REMARKS 

nowfieid-Dependent  Mixed  ExpUcit- 
w  ^  method  lias  be«i  inmuiacsd  for  both  aooieacting  and  rsaoting 
Tils  method  is  b^eved  to  be  panicnlarly  nseftd  for  snpersonic  combnstiom 

““  ““  “““  implicitness  pa^eters  and  var^^ 
eJmtions  of  Damimhlet  nnmber  as  relaud  to  the  source  term  Jacobians  play  a 
^ecant  role  in  the  hnite  rate  chemistry,  especiaUy  with  shock  wa,eTurb,!!ent 
^no^  layer  interactions  Chemical  reactions  are  either  enhanced  or  dimisnished 
pnaence  of  comply  Sow  physics  including  transition  from  lamiei.r  to  turbul- 
^  from  incompressible  to  compressible,  and  from  invisdd  to  viscous  flows,  mid 
resniting  m  the  vanous  degree  of  efficiency  in  combustion. 

■  ■  ■  ettehsive  research  wiU  be  required  before  the  full  assessment  of  tie 

TliUna  MEI  slieme  for  applications  to  supersonic  combustion  can  be  made. 

tote  T’  o'  preliminary  results,  eitpected  certaiiy 

to  be  reimed  and  reenforced  m  the  future  research.  ^ 
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9)  h) 

Fifure  7.  Various  flowfield  contours  for  transverse  hydrogen  fuel  injection,  (a)  Density 
contours  (frozen);  max  =  0.6,  min  =  0.058,  A  =  0.021  kg/m^.  (6)  Density  contours 
(reacting):  max  =  1.813,  min  =  0.0842,  A  =  0.049  kg/m^.  (c)  Pressure  contours  (reacting): 
max  =  0.4053,  min  =  0.0017,  A  —  0.025  MPa.  (d)  Temperature  contours  (reacting): 
max  =  2982,  min  =  345,  A  ss  195  K.  (e)  Hj  mass  fraction:  max  =  1.0,  min  =  0.0, 
A  ss  0.024.  (/)  Oj  mass  fraction:  max  =  0.095,  min  =  0.0,  A  =  0.002.  [g)  OH  mass 
fraction:  max  =  0.0602,  min  =  OJ),  A  s  0.0024.  (A)  H2O  mass  fraction:  max  =  0.0081, 
min  s  0.0,  A  =  0.0003 
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Abstract 


Despite  significant  achievements  in  computational  fluid  dynamics,  there  still  remain  many  fluid  flow  phenomena  not  well  understood.  For 
example,  ihe  prediction  ot  tempeniture  disinbuiions  is  inaccurate  when  temperature  gradient.';  are  nigh,  particularly  in  shock  wa\  e  lurbuleni 
boundary  layer  interactions  close  to  the  wall.  Complexities  oi  fluid  flow  phenomena  include  transition  to  turbulence,  reiaminanzatton. 
separated  flows,  transition  between  viscous  and  inviscid.  incompressible  and  compressible  flows,  among  others,  in  ail  speed  regimes.  The 
purpose  01  this  paper  is  to  introduce  a  new  approach,  called  the  Flowjield-Dependeni  Mixed  Explicit-Implicit  i  FDMEI)  metnod.  in  an 
attempt  to  rtsolve  these  difficult  i.ssues  in  CFD.  In  this  process,  a  total  of  six  implicitness  parameters  characteristic  of  the  current  flowheld 
are  introduced.  They  are  calculated  from  the  current  flowheld  or  changes  of  Mach  numbers.  Reynolds  numbers.  Peclet  numbers,  and 
amkohler  numbers  (if  reacting)  at  each  nodal  point  and  time  step.  This  implies  that  every  nodal  point  or  element  is  provided  with  different 
or  unique  numencal  scheme  according  to  their  current  flowheld  situations,  whether  compressible,  incompressible,  viscous,  inviscid.  laminar, 
turbulent,  reacting,  or  nonreacting.  In  this  procedure,  discontinutties  or  fluctuations  of  all  vanables  between  adjacent  nodal  points  are 
determined  accurately.  If  these  implicitne.ss  parameters  are  fixed  to  certain  numbers  instead  of  being  calculated  from  the  fiowfield 
information,  then  practically  all  currently  available  schemes  of  finite  differences  or  finite  elements  anse  as’  special  cases.  Some  benchmark 
problems  to  be  presented  in  this  paper  will  show  the  validity,  accuracy,  and  efficiency  of  the  proposed  methodology. 


1.  Introduction 


Nearly  half  a  century  has  elapsed  since  the  digital  computer  revolutionized  computational  technologies  in 
engineering  and  mathematical  physics.  During  this  time  finite  difference  methods  iFDM)  have  dominated  the 

dynamics  (CFD)  [1-7],  whereas  the  opposite  is  true  for  finite  element  methods 
(FEM)  in  solid  mechanics.  In  recent  years,  however,  the  trend  toward  finite  element  methods  in  CFD  appears  to 
be  increasingly  favorable  (8-14). 

In  general,  the  analyst  preoccupied  with  the  methods  of  his  choice  based  on  his  educational  background  or 
research  experience  is  seldom  motivated  to  investigate  other  options.  Thus,  today  the  gap  between  these  two 
«mp  ines  is  widely  apart,  despite  the  fact  that  the  thorough  understanding  of  the  relations  between  FDM  and 
miM  purpose  of  this  paper  is  an  attempt  to  call  for  a  new  approach  in  which  both  FEM  and 

TOM  can  be  united  toward  the  common  goal  of  achieving  the  highest  level  of  accuracy  and  efficiency  in  CFD. 
Similarities  and  dissimilarities  must  be  identified  in  order  to  recognize  merits  and  demerits  of  each  method  and 
to  enable  the  analysts  to  choose  the  most  desirable  approach  suitable  for  the  particular  task  at  hand. 
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One  ot  the  most  imponant  questions  m  CFD  i>  now  to  aeai  u  un  iarge  graments  or  the  vanaoie  idensitv, 
\eocit\.  presst^e.  temperature,  and  source  termsi.  Raoid  chances  of  Mach  numoers.  Revnoids  numoers  Peclet 
nuinbers.  and  Damkohler  numbers  t.f  reacting)  oetween  adjacent  nodal  points  or  elements  can  be  a 'crucial 
actor  m  determining  whether  the  chosen  computational  scheme  will  succeed  or  fail.  Funhermore  Drooer 
treatments  tor  incompressibility  and  compressibility,  viscous  and  inviscid  flows,  subsonic  and  supersonic  flows 
laminar  and  turoulent  flows,  nonreacting  and  reacting  flows  are  extremely  imponant.  The  most  eeneral  case  of 
fluid  dynamics  where  these  various  flow  propenies  may  be  depicted  m  external  and  internal  hypersonic  flows  is 
shown  m  Fig.  l(a.b).  A  typical  reacting  flow  t hydrogen-air  reaction)  can  also  be  seen  in  Fig  1(c) 

computer  code  be  made  available  to  satisfy  all  the  requirements  mentioned 
hr.fh  ^^ihematicai  formulation  lead  to  most  of  the  cunently  available  computational  schemes 

j  ^  “  special  cases  .’  Most  imponantly.  will  such  an  approach  guarantee  accuracy  and 

exmpirproblems  qiiestions  positively,  based  on  the  results  obtained  through 

Toward  this  go^.  our  approach  is  based  on  the  following  procedure  [15.16],  known  as  the  Flowtield- 
Dependent  Mixed  Explicit-Implicit  (FDMEI)  scheme:  riowneia 

(a)  Write  the  Navier-Stokes  system  of  equations  in  a  conservation  form. 


\.r.  Yoon  t'l  ill.  C:innui.  .\ff!nna\  ^nn;  r.r-jr?  ..'o!  fyX)-()()() 

tb)  Expana  the  consen'aiion  vanaole  C  n  Tavior  senes  ud  to  and  inaiudinc  tne  secona-oraer  time 
derivaiives  or  the  conservation  variables. 

(c)  Introduce  in  step  (b)  six  different  tiowrieid-depenaent  implicitness  parameters  which  are  calculated  from 
the  changes  m  Mach  numbers.  Reynolds  numoers.  Peclet  numoers.  and  Damkohler  numbers  (if  reacuno 
between  nodal  points  or  local  elements. 

(d)  Substitute  step  tai  and  (c)  into  step  tb)  to  obtain  the  increments  of  the  conser\'ation  variables  AL'""’  .\s  a 
result,  the  final  form  resembles  the  implicit  factored  scheme  of  Beam  and  Warming  [1],  but  much  more 
rigorous. 

(e)  Step  (d)  may  be  used  either  in  FDM  or  FEM. 

The  computational  procedure  as  described  above  is  capable  of  resolving  complex  propenies  of  fluid  flows  in 
general  with  shock  waves,  mrbulence.  and  reeling  flows  in  particular. 

(1)  Shock  waves  in  compressible  flows  are  dependent  on  changes  in  Mach  number  between  nodal  points  in 
FDM  and  within  local  elements  in  FEM.  Shock  wave  discontinuities  are  characterized  by  these  changes 
in  Mach  number. 

(2)  Incompressible  turbulent  flows  are  dependent  on  changes  in  Reynolds  number  between  nodal  points  in 
FDM  and  within  local  elements  in  FEM.  Incompressibility  conditions  are  characterized  by  these  changes 
in  Reynolds  number. 

(3)  Compressible  turbulent  flows  are  dependent  on  changes  in  both  Mach  number  and  Reynolds  number 
between  nodal  points  in  FDM  and  within  local  elements  in  FEM.  Dilatational  dissipation  is  characterized 
by  these  changes  in  Mach  number  and  Reynolds  number. 

(4)  Hi§h  temperature  sradient  flows  are  dependent  on  changes  in  Peclet  number  between  nodal  points  in 
FDM  and  within  local  elements  in  FEM.  The  convection  vs.  diffusion  in  heat  transfer  is  characienzed  by 
these  changes  in  Peclet  number. 

(5)  Reaction  flows  are  dependent  on  changes  in  Damkohler  number  between  nodal  points  in  FDM  and  within 
local  elements  in  FEM.  The  mass  source  vs.  convective  transfer,  mass  source  vs.  diffusive  transfer,  heat 
source  vs.  convective  heat  transfer,  heat  source  vs.  conductive  heat  transfer,  and  heat  source  vs.  diffusive 
heat  transfer  are  characterized  by  these  changes  in  Damkohler  number. 

(6)  Direct  numerical  simulation  (DNS)  in  which  mesh  refinements  are  carried  out  until  turbulence  length 
microscales  are  resolved  without  turbulence  models  can  not  be  reliable  particularly  for  high  speed 
compressible  turbulent  flows  unless  the  computational  .scheme  is  capable  of  treating  high  gradients  of 
variables  as  described  in  (3)  above.  To  improve  turbulence  calculations,  Legendre  polynomial  spectral 
modes  may  be  added  as  shown  in  [15].  Whether  or  not  the  spectral  mode  approach  is  advantageous  for  an 
overall  computational  efficiency  remains  to  be  seen.  Due  to  the  limitation  of  computer  time,  the  example 
problems  in  this  paper  are  not  intended  for  DNS  micro.scaie  resolutions. 

Details  of  the  mathematical  formulations  as  described  above  are  presented  in  Section  2,  implementation  and 
computational  process  m  Section  3,  some  example  problems  in  Section  4.  and  concluding  remarks  in  Section  5. 


2.  Mathematical  formulations 

For  the  general  purpose  program  considering  the  compressible  viscous  reacting  flows,  we  write  the 
conservation  form  of  the  Navier-Stokes  system  of  equations  as 

dU  DF.  dG; 

(it  (lx I  dx I  ^  (n 

where  £/,  F.,  G,  and  B  denote  the  conservation  flow  variables,  convection  flux  variables,  diffusion  flux  variables, 
and  source  terms,  respectively. 
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wnere  r  -  tne  cnemicai  species.  H''  is  the  zero-point  entnaipy.  u  .  is  the  reaction  rate,  ana  D, :> 

the  binary  ditfusivity.  Additional  eauations  for  vibrational  and  electronic  energies  mav  be  included  mil)  tor 
hypersonics. 

Expanding  the  conserv'ation  variables  U  in  Tavlor  senes  includinc  the  rirst  and  second  derivatives,  we  have 


U"^'  =U"  ^  \t—. - —  - 


-O(At') 


where  and  j,  are  the  implicitness  parameters  defined  such  that 
dU"^''  dU"  dAU"^' 


d-u"""-  d'y"  d-AU" 
dY  ~  dr  dt 


0«s,  ^  1 


with  AU  -U  *  '  -  {/“,  It  is  assumed  that  the  convection  flux  is  a  function  of  U  and  the  diffusion  flux  G 
is  a  function  of  both  U  and  its  gradient  (7^.  Thus,  we  have 

m  __  rlF,  aG 

(5) 


d'U  _d_ 

dr  '^•v 


d  (  hU\  h  (  dU\  d-  (  dU\  /  , 
LV;  v“'  dt  )  dx.  V‘  dt  )  dx.  dXj  V"  dt  / 


where  the  convection  Jacobian  a,,  the  diffusion  Jacobian  the  diffusion  gradient  Jacobian  c,^,  and  the  source 
Jacobian  d  are  defined  as  ‘ 

dF/  dG;  dG,  dB 

‘^'^dU'  ^'~dU' 

Substituting  (3)-(6)  into  (2)  and  assuming  the  product  of  the  diffusion  gradient  Jacobian  with  third-order  spatial 
derivatives  to  be  negligible,  we  obtain 


dAF""'  dAG""' 


2  d, 


dF"  dG'/ 

dx-  dx 


\  (dF"  dG" 

-fl"  - —-B" 


dx  dx 


dAFT'  dAG"*' 


dAFr'  dAG" 


-AB"*' 


+  0(Ar') 


where  all  Jacobians  are  considered  to  be  constant  within  an  incremental  numerical  time  step,  but  allowed  to  be 
updated  at  consecutive  time  steps. 

In  order  to  provide  different  implicitness  (different  numerical  treatments  or  schemes)  to  different  physical 
quantities,  we  reassign  5 ,  and  s,  associated  with  the  diffusion  and  source  term,  respectively, 

s ,  AG,  => Xj  AG, .  s,  AB  =»s,  AB  (9a) 

Sj  AG,  =>  AG, .  s.  AB  AB  (9b) 


with  the  various  implicitness  parameters  defined  as 
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y,  =  rirst-oraer  convection  impiicuness  parameter 
=  second-order  convection  impiicuness  parameter 
5,  =  tirst-order  diffusion  implicitness  parameter 


=  second-order  diffusion  implicitness  parameter 
=  rirsi-order  source  term  implicitness  parameter 

Sf^  =  second-order  source  term  implicitness  parameter 

The  first  order  implicitness  parameters and  will  be  shown  to  be  fiowrieid  dependent  with  the  solution 
accuracy  assured  by  taking  into  account  the  flowfieid  gradients,  whereas  the  second  order  implicitness 
parameters  and  which  are  also  flowfieid  dependent,  mainly  act  as  artificial  viscosity,  contributins  to  the 
solution  stability. 

Substituting  these  implicitness  parameters  as  defined  in  (9)  into  (8).  we  obtain 


A/-  /  \ ^hafl.■^ba,)^U 
2 


I 


db,AU"'' 


dr 


-d 


AO:  AU‘ 

OX; 


dXi  dx- 


-S;d  A£/""' 

d\a:bf+bb.)AU‘'^' 


dX;  dx- 


(  db  ^U"*'  d'c..^u’" 

-d{- 


i  \ 


\ 


dx  dx 


with 


AB 


ir 

1 


\ 


dia.  -H  b  )  AL^"' 
fix 

fiF"  //G" 


4 


f)F'‘ 

iix 


■,G" 


‘I  ^  -  O"  \  /  .^F"  :,G"  \ 

®  r  "I  ij- - — - ®“)  J  - 


dx 


0(A/')  =0 


dB 

dU 


At/"^'  =dAU' 


Rearranging  (10),  we  obtain 
'  dX: 


R=A^U‘' 


(£,  At/"’ ' )  +  (E,  At/" ’ ' )  +  Q  '  +  0(A;' 


or 


( 


with 


Ar" 


A  ^  I  H-  Ar  s^d  d — —  s^d 


A/ 


E,  =  +  Sybi)  +  —  + 1», )  +  5,  da,  +  s,db,] 

_  ,  Ar" 

E,}  =  Ar  s^c,^  -  —  (jj(a,a.  +  ba.)  +  s,(ap^  +  b,hj  -  dc,.)\ 


(10) 

(11) 

(12) 

(13) 

(14) 

(15) 

(16) 


(17) 


An  alternative  scheme  is  to  allow  the  source  term  in  the  LHS  of  ( 13)  to  lag  from  «  -r  1  to  n  so  that  (13)  may  be 
written  as 
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w,/  .  UJO-7)  ()0()^(y^) 


-£  -■ 


rr  ax 


£  ,  )  A£" 


=  -i)" 


;  18) 


with 


^  ^^d)(F':  -  G'r,  ^  a,  ~b,  ,B" 

El  =  111  A/(a,  +  Z>, ) .  with  m  =  e/([  +  A) 

(20) 


Eij  =  m^tc.j 


2"  = 


At 


/dF" 


bG'!\  ^ 

~)+~AU-' 


(21) 


(22) 


1  +  f  \  ajC; 

d^ZluTc  ZZl  ort^rLHs"™hrBTrWa™.nG  ‘""‘“dec!  .n  the  seconc 

.f  we  set  =  ,n.  ,,  ",  =  .  IT-ZT1  Tf  equat.on.< 

analysts  of  the  ^ ^ 

parameter //;  to  be  0.639  ^ //I  0  75  It  rnn  eh  '  u  ‘  ^  ^  3  wiii  fix  the  implicitness 

capable  of  producing  practically  all  existine  FDM  FDMEI  equations  as  derived  in  (13)  or  (18)  are 

A.  ^  ^  PDM  and  FEM  schemes.  Some  examples  are  shown  tn  Appendix 

the  current  flowfield  variable™!  each  ^Tevenr  ^daf^  approach  is  to  obtain  the  implicitness  parameters  from 
certain  predetermined  numbers  and  using  them  f  fixing  the  implicitness  parameters  to 

vanation  from  one  point  to  another  S  iZrZ 

bases  as  depicted  in  Fig.  2  The  final  values  of  '  parameters  may  be  determined  for  spatial  and  temporal 
obtained  as  the  average  of  both  spatial  and  tempo'll  ertriburions'"''"''  """ 

Convection  implicitness  parameters 


■s,  = 


minir,  1) 
0 


r>  a 

M..„=0 


A  ^  —  s, 


0<n  < 


with 


'■=VM;„-m^..„/m„ 


(23) 


(24) 


wlSb  rLTfi^relemenTm  S  hZLZaTiZZt  P°‘"“  ° 

»  »d  »  -  I  for  impUctos/^^Cr^T r,™  T7'“”  “ 

Here,  It  IS  seen  that  s,  is  directly  related  to  ih^  n  ^  «  !!!  ’  ^  user-spectfied  small  number  a  =0.01) 

primary  role  of  a,  ,r  m  erniumtS  “  J  «"  “■="  'fx  «=->;.  Th, 

whereas  that  of  s.  is  ,0  aa  as  anifidai  viscosity,  to  mlutimSlit“°"'"’°‘^“"®  Sratiienis 

Diffusion  implicitness  parameters 


s,  =i 


minis.  1 )  s>  fi 

0  s<P,  Re,„,„  =  0 .  or  Pe 

‘  RCm..  =0. 


7=0 


or  Peri...  =0 


S4  =  s',' .  0  <  /I  <  1 


(25) 
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Fig.  2.  Spatial  and  temporal  flowfteid  dependent  implicitnes.s  parameienc  (a)  Idealized  turbulence  length  scales  assumed  to  be  wuhm  each 
elemeni-spatialiy  evolving  turbulence.  The  maximum  and  minimum  value  of  M.  Re.  Pc  and  Daare  those  among  the  comer  nodes  within  an 
clement,  (b)  idealized  turbulence  time  scales  a.s.sumed  to  be  within  each  time  step-temporal ly  evolving  turbulence.  The  maximum  and 
minimum  value  of  M.  Re,  Pc  and  Da  are  identified  at  time  steps  n  and  /i-l. 


with 

^=VRe;„  or  J  =VPeL  "  Pe^i„/Pe„,„  (26a.b) 

where  the  maximum  and  minimum  Reynolds  numbers  or  maximum  and  minimum  Peclet  numbers  are  calculated 
simiiarily  as  in  5,  for  spatial  and  temporal  implicitness  parameters,  and  /3  is  a  user-specified  small  number 
iP  “0.01).  If  temperature  gradients  are  large,  it  is  possible  that  Peclet  numbers  instead  of  Reynolds  numbers 
will  dictate  the  diffusion  implicitness  parameters.  The  larger  value  of  .v,  is  to  be  chosen,  as  obtained  either  from 
(26a)  or  i26b).  Note  also  that  5^  =  with  ensuring  the  .solution  accuracy  by  taking  into  account  the  diffusion 
gradients,  and  here  again,  5^  plays  the  role  of  artificial  viscosity,  for  solution  stability. 

Source  term  implicitness  parameters 

For  the  case  of  chemically  reacting  flows  the  Dq 

'min(/,  1 )  t^y 

U  Da„i„=0 

with 

t  =VDa;„  -  Da^„/Da^„  (28) 

where  the  maximum  and  minimum  Damkohier  numbers  are  calculated  similarly  as  in  5,  and  for  spatial  and 
temporal  implicitness  parameters,  and  y  is  a  user-specified  small  number  (y  =  0,01).  The  relationship  between 
55  and  is  similar  to  those  for  convection  and  diffusion  implicitness  parameters  such  that  =  55  with  and  S(^ 
controlling  the  solution  accuracy  and  solution  stability,  respectively’.  The  average  of  both  spatial  and  temporal 
implicitness  parameters  will  be  adopted  for  use  in  computations  at  any  point  (element)  and  time. 


(Damkohier  number)  must  be  used 

0  <  rt  <  I  (27) 


A  T".  Ynoii  et 


ntwut.  Mt'tnnas  Anni.  Mccn.  £/u^rp.  tJni  i  1997)  OOO-^f'HX) 


between  ail  physical  phenomena  and  the  correspondins  numencai  treatments  are  cnaractenzea 

parameter^rV""  second-orXT^rness 

to  orovirip  nb'  '  computational  accuracy  and  computational  stabilitv,  respectiveiv  The  idea  is 

coi^putation^Turv  “L^theVet  ’  T"' 

modified  from  those  reponed  in  115  161  m  r'  ^econd^order  implicitness  parameters  have  been 

thoueht  that  ,hP  c!  I  .  ^  requirements  (Fia.  3).  Initially  it  was 

first-mder  implicitness"  7^L!^fT=l-7T-l 

ihTfi^sST^^^^^^^  -pec^;iy:fo?the  aid 

for  the  high  values  ofThe'&nrteTimrtvr”'’''  "“hierical  viscosities 

secoitd^der  impltchtess  parameters  Tern  pmvTdeTTchTa.  '  fo'cT”"'  ''“““  "" 

[15],  However,  it  was  noted  that  hnrh  fir«  ^  -  r ,.  0.5).  etc.  as  expenmented  in 

both  extremes  at  aero  and  unim  »hh  S  se»,^ordt°'^ 

values  of  rhe  hrst-nme,  lmrii...L.  -.-I.  pemmeieia  being  reasonably  large  for  all 

the  nonlinear  conttnuous  functions  onTfot  o^  ■rntfic'!'  ■”Plici.ness  parameter  given  above  are 

range  of  the  constant  n  is  0</i  <  l  althn  i,  |  parameters  sansfying  these  requirements.  The 

convergence  rate  for  reasonably  high  CFI  h  '  optimum,  exhibiting  the  best 

The  fiowfield  denendent  imnli  numbers  in  the  example  problems  presented  in  Section  4. 

schemes  ,„  be  automabcally  general  rs^^Sd^  i“ 

waves  and'^tu^iiT*^*''^*''!'^''  •'i  •^3  control  all  high  gradient  phenomena  such  as  shock 

Rcyndl  (m  Pecren  ■  T  of  local  Mach  numSr^  anS 

flowfields  The  contoureof  indicative  of  the  actual  local  element 

anH  c  h»  I  /  ?  ^  ^  parameters  closely  resemble  the  flowfields  themselves  with  both  j 

equate  computational  siability  (anificial  viscosity)  as  they  were  originally  introduced 


A.  .ig -1-1, .,,-1-,^ 

a  *3  -  Bwtfl  -  g,,  0.5). -  na*<i  -  gj,  OJ).  g,  -  mnt I  -  g,.  0.5) 

^  =r;.f,  =g,*.0<a<t 


optimum  at =  i  for  the  secLd^rdel  impli^ess  plramew^'m"''*  ’'“‘“‘'“''s  occur  in  the  range.  0  <  /i  <  I.  with  i 

panimeten.  *’  *°  preserve  the  solution  accuracy  as  dictated  by  the  first-order  tmpiicitne: 
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mio  the  secona-oraer  time  aenvative  term  or  me  Tayior  senes  expansion  or  the  conservation  riow 
vanabies  U  The  primarx'  roie  of  .v.  and  s.  to  provide  comouiaiionai  stabiliiv. 

f  3)  The  s,  terms  represent  convection.  Tliis  implies  mat  \t  s,  =  0  then  tne  effect  or  convection  is  smaii.  The 
compuiationai  scheme  is  automatically  altered  to  take  this  effect  into  account,  with  the  govemma 
equations  being  predominantly  parabolic-elliptic.  Note  that  these  effects  are  confined  at  not  at  U" . 

(4)  The  terms  are  associated  with  diffusion.  Thus,  with  =  0.  the  effect  of  viscosity  or  diffusion  is  small 
and  the  computational  scheme  is  automatically  switched  to  that  of  Euler  equations  where  the  govemma 
equations  are  predominantly  hyperbolic. 

(5)  If  the  hrst-order  implicitness  parameters  and  are  nonzero,  this  indicates  a  typical  situation  for  the 

mixed  hyperbolic,  parabolic  and  elliptic  nature  of  the  Navier-Stokes  system  of  equations,  with 
convection  and  diffusion  being  equally  important.  This  is  the  case  for  incompressible  flows  at  low 
speeds.  The  unique  property  of  the  FDMEI  scheme  is  its  capability  to  control  pressure  oscillations 
adequately  without  resorting  to  the  separate  hyperbolic  elliptic  pressure  equation  for  pressure 
corrections.  The  capability  of  FDMEI  scheme  to  handle  incompressible  flows  is  achieved  by  a  delicate 
balance  between  r,  and  as  determined  by  the  local  Mach  numbers  and  Reynolds  (or  Peclet)  numbers. 
If  the  flow  is  completely  incompressible  (M  =  0),  the  criteria  given  by  (23)  leads  to  5,  =  1.  whereas  the 
implicitness  parameter  Sy  is  to  be  determined  according  to  the  criteria  given  in  (25).  Make  a  note  of  the 
presence  of  the  convection-diffusion  interaction  terms  given  by  the  product  of  b  a.  in  the  s.  terms  and 
apj  in  the  terms.  These  terms  allow  interactions  between  convection  and  diffusion  in  the  viscous 
incompressible  and/or  viscous  compressible  flows. 

(6)  If  temperature  gradients  rather  than  velocity  gradients  dominate  the  flowheld,  then  is  governed  by  the 
Peclet  number  rather  than  by  the  Reynolds  number.  Such  cases  arise  in  high  speed,  high  temperature 
compressible  flows  close  to  the  wail. 

(7)  In  the  case  of  reacting  flows  the  source  terms  B  -contain  the  reaction  rates  which  are  functions  of  the 
flowfieid  variables.  With  widely  disparate  time  and  length  scales  involved  in  the  fast  and  slow  chemical 
reaction  rates  of  various  chemical  species  as  characterized  by  Damkohier  numbers,  the  first-order  source 
term  implicitness  parameter  is  instrumental  in  dealing  with  the  stiffness  of  the  resulting  equations  to 
obtain  convergence  to  accurate  solutions.  On  the  other  hand,  the  second-order  source  term  implicitness 
parameter  contribute  to  the  stability  of  solutions.  It  is  seen  that  the  criteria  given  by  (27)  will  adjust 
the  reaction  rate  terms  in  accordance  with  the  ratio  of  the  diffusion  time  to  the  reaction  time  in  finite  rate 
chemistry  so  as  to  assure  the  accurate  solutions  with  computational  stability. 

(8)  Various  definitions  of  Peclet  number  and  Damkohier  numbers  (Table  1)  between  the  energy  and  species 
equations  should  be  checked.  Whichever  definition  provides  larger  values  of  Sy  and  must  be  used.  The 
summary  of  the  above  definitions  for  implicitness  parameters  is  shown  in  Table  2. 

(9)  The  transition  to  turbulence  is  a  natural  flow  process  as  the  Reynolds  number  increases,  causing  the 
gradients  of  any  or  all  flow  variables  to  increase.  This  phenomenon  is  the  physical  instability  and  is 
detected  by  the  increase  of  Sy  if  the  flow  is  incompressible,  but  by  both  Sy  and  5,  if  the  flow  is 
compressible.  Such  physical  instability  is  likely  to  trigger  the  numerical  instability,  but  will  be  countered 
by  the  second-order  implicitness  parameters  and/or  to  ensure  numerical  stability  automatically.  In 
this  process,  these  flowfieid  dependent  implicitness  parameters  are  capable  of  capturing  relaminarization. 
compressibility  effect  or  dilatational  turbulent  energy  dissipation,  and  turbulent  unsteady  fluctuations. 

(10)  An  important  contribution  of  the  first-order  implicitness  parameters  is  the  fact  that  they  can  be  used  as 
error  indicators  for  adaptive  mesh  generations.  That  is,  the  larger  the  implicitness  parameters  the  higher 
the  gradients  of  any  flow  variables.  Whichever  governs  (largest  first-order  implicitness  parameters)  will 
indicate  the  need  for  mesh  refinements.  In  this  case,  all  variables  (density,  velocity,  pressure, 
temperature,  species  mass  fraction)  participate  in  resolving  the  adaptive  mesh,  contrary  to  the 
conventional  definitions  of  the  error  indicators. 

(11)  Physically,  the  implicitness  parameters  will  influence  the  magnitudes  of  Jacobians.  Thus,  Item  8  above 

may  be  modified  so  that  the  diffusion  implicitness  parameters  Sy  and  as  calculated  from  Reynolds 
number  and  Peclet  number  can  be  applied  to  the  Jacobians  corresponding  to  the  momentum 

equations  and  energy  equation,  respectively.  Furthermore,  two  different  definitions  of  Peclet  number 
(PCp  Pe,j)  would  require  the  ^3  and  5^  as  calculated  from  the  energy  and  species  equations  to  be  applied 
to  the  corresponding  terms  of  the  Jacobians.  Similar  applications  for  the  source  term  implicitness 
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Table  t 

Dertnmons  or  nonaimcnsionai  tiowneld  auunmies 


•  V  If  -  -\n  ^  U  ^  T  V  -  r ij 

A  B  ^ 

v-pvj^  c,_  jr-T-  ;,o  VKj^  t^dr-v-jvr  = 

y 


V ■  ( py] V )  -  V  ( po  vy. }  u- 


N 


Mach  number 

M 

u 

A  incnial  force 

a 

B  pressure  force 

Reynolds  number 

Re 

puL 

A  initial  force 

C  viscous  force 

Peclci  number.  I 

Pc, 

E  convective  heat  transfer 

k 

G  conductive  heat  transier 

Peclei  number.  [1 

Pe„ 

uL 

I  convective  mass  transfer 

w,, 

"d 

i  diffusive  mass  iranster 

Dumkohler  numocr.  I 

Da, 

Lw. 

K  mass  source 

puY, 

I  convective  itranster 

Damkohler  number.  11 

Da„ 

K  mass  source 

pDY, 

J  diffusive  (transfer 

Damkbhier  number.  Ill 

Da„, 

N  heat  source 

Hu 

E  convective  heat  iranster 

Damkohler  number.  IV 

Da,v 

sE. 

H  heat  source 

kT 

C  conductive  heat  transfer 

Damkbhier  number.  V 

Dav 

SiL 

N  heat  source 

HD 

F  diffusive  heat  transfer 

parameters  s,  and  should  be  followed  for  the  source  term  Jacobian  d.  In  this  way,  high  temperature 
gradients  arising  from  the  momentum  and  energy  equations  and  the  finite  rate  chemistry  governed  by  the 
energy  and  species  equations  can  be  resolved  accordinslv. 

Ihe  t-DMEl  equations  as  given  in  il3)^may  be  solved  by  either  FDM  or  FEM.  The  standard  linear  Galerkin 
Wroximations  ot  FEM  lead  to  the  results  of  central  differences  of  FDM.  However,  the  main  difference  between 
TOM  and  FEM  anses  when  integration  by  parts  is  performed  in  FEM  and  the  explicit  terms  of  Neumann 
toundary  conditions  naturally'  appear  as  boundary  integral  forms.  Thus,  all  Neumann  boundary  conditions  can 
irecily  specified  at  boundaries  in  FEM.  This  is  not  the  case  for  FDM.  Often,  a  rather  cumbersome  process 
must  be  taken  for  Neumann  boundary  conditions  in  FDM. 

When  dealing  with  all  speed  flow  regimes  such  as  in  shock  wave  turbulent  boundary  layer  interactions  where 
compressible  and  incompressible  flows,  viscous  and  inviscid  flows,  and  laminar  and  turbulent  flows  are 
intermingled  throughout  the  flowfield  domain,  a  computational  scheme  intended  for  only  one  type  of  flow 
p  ysics  and  that  does  not  account  for  other  types  of  flow  phenomena  will  fail.  For  example,  the  flow  close  to  the 
wall  in  shock  wave  turbulent  boundary  layer  interactions  is  incompressible  (M^O.l).  whereas  away  from  the 
wall  the  flow  is  compressible  (supersonic  or  hypersonic).  In  this  case,  viscous  flows  change  to  inviscid  flows.  In 
between  these  two  extremes  the  flowfield  changes  continuously,  oscillating  back  and  fonh  across  the  boundary 
layers  of  velocity  and  entropy,  and  leading  edge  and  bow  shocks.  At  any  given  computational  nodal  point  or 
element,  padients  of  each  variable  (density,  pressure,  velocity  and  temperature)  may  be  very  small  or  very 
arge,  so  large  that  practically  all  currently  available  computational  methods  may  fail.  In  order  to  succeed,  it  is 
n^essary  that  the  current  flow  physics  everywhere  be  identified  and  so  recognized,  with  specific  computational 

sc  ernes  accorded  to  each  and  every  computational  nodal  point  and  element.  It  is  clear  that  such  accommoda¬ 
tions  are  available  in  (13)  or  (if). 
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Table  2 

Rowtieid  deoenaent  imDiicimess  oarameters 

Convection  gradient  behavior  - - 

'  — Rr.Tf-orrtt'r  convection  impitcitness  parameter 

Ensures  solution  accuracy 

Strongly  fiowheld  dependent,  witn  hign  graoieni.s  cnar- 
actenzed  by  large  Mach  numocr  cnanees  ociween  nodal 
points  or  within  eiemeni  and  between  time  steps 


{mint/-,  I  )  r  >  a 
0  r<a. 

*  =0 

~  Vf  /M 


5 , ^Second-order  convection  implicitness  parameter 


S:=s'l.  0<«<l 


Diffusion  gradient  behavior. 

'  ~^Eirsf~ot  (ier  diffusion  implicitness  parameter 


Ensures  solution  stability 

Rowfield  dependent  artificial  viscosity  tor  convection 
process 


{mints.  1 1  V  >  ^ 

0  s</3,  Re,„,„  =«0,  orPe,„,^?^0 

*  or  Pe„„,  =  0 

'  =VReL  -  Re:„,./Re,„,.  or  t  =\/Pe;„  -  Pei^/Pe,,,,, 


s  ^--Second-order  diffusion  implicitness  parameter 

54  0</i<| 

Source  term  gradient  behavior 
s ,—Firsr-oraer  source  term  implicitness  parameter 


Ensures  solution  accuracy 

Strongly  flowficld  dependent,  with  high  gradients  char¬ 
acterized  by  large  changes  in  Reynolds  number  or  Peciet 
number  between  nodal  points  or  within  element  and 
between  time  steps.  Diffusion  gradient  behavior  may  be 
dictated  by  Peciet  number  when  temperature  gradients  arc 
high.  Choose  which^cver  <  Reynolds  or  Peciet  number! 
provides  the  larger  value  for  s, 


Ensures  solution  stability 

Flowheld  dependent  artificial  viscosity  for  diffusion  pro¬ 
cess 


{min(r.  n  t  ^  a 

0  r<a.Da *0 

'  Da,.,.=0 

-•  =  VOai,  -  Dai,. /Da . 


Ensures  solution  accuracy 

Strongly  fiowfield  dependent,  with  high  reaction  rate 
gradients  choractenzed  by  large  Damkbhler  number 
changes  between  nodal  points  or  within  element  and 
between  time  steps 


s ^Second-order  source  term  implicitness  parameter 


5^  =  j;' .  0  <  «  <  I 


Ensures  solution  stability 

Rowfield  dependent  anificial  viscosity  for  reaction  process 


3.  Implementation  and  computational  process 

governing  equations  for  the  Taylor  series-modified  Navier— Stokes  system  of  equations, 
or  to  the  finite  volume  method  (FVM).  For  FDM  applications  the 
or  er  an  secon  order  spatial  derivatives  may  be  written  in  central  difference  schemes.  For  example,  we 
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may  write  tor  any  variable  ii 
-r.,  ^ 

rt.rl,.,  ~ 


d'u 

dx'- 


Z\x 

*  I  ,  -  2u  .  ^ 

I? 


29) 


(30) 


with  analogous  formulas  for  the  ,v  derivatives  in  two  dimensions  or  corresponding  extensions  to  three- 

dimensional  problems.  Any  other  difference  schemes  of  higher-order  accuracy  mav  also  be  chosen  as  deemed 
necessary.  '  ' 

For  applications  to  FEM  we  begin  by  expressing  the  conservation  and  flux  vanables  and  source  terms  as  a 
linear  combination  ot  trial  functions  <P^  with  the  nodal  values  of  these  variables. 

U(x,  t)  =  0„(x)f/„(/) ,  F,(x,  r)  =  0„(r)F„,(f) 

G,(x,  t)  =  0„(x)G„,(O  .  B{x,  t)  =  0„(x)B„(r) 

Applying  the  generalized  Galerkin  approximations  to  (11)  we  obtain 
<P.RiU.F  .G,.B)An=0 

Jn 


or 


ag;:'  =  w;,  +  n':^ 


where 


(31) 


(32) 


(33) 


%  d/2  .  77„  -  5,,  -  ( A/  Sj  +  -^  (34) 

"  fa  +  hK,)  +  ^  +  s,dja„^  +  hj  + 

-  —  U3(a,>,u^„  +  -I-  sMa,b„.  ^  ^  d/2 

+  -  —  [sj(a,„a,„  +  +  b,„b„^  -  J„.  dF  (35) 

“  fa  ^a->)  +  — -h  G"j+Y{a,,,  + 

-  —  (a,>,  +  b,„)(F;j,  +  G"^..)0„ ,  <Pg^  -H  [a/  f;,  +  d/2  06) 

=  fr  ^  .  G;.J  -  ^  (a.,.  ^  b,„)B;]  0.0, 

+  -T  +  g;.J0.0,  J  dF  (37) 

♦ 

where  0.  represents  the  Neumann  boundary  tr«l  and  test  functions,  with  a.  /3  denoting  the  global  node  number 
the  number  of  conservation  variables  at  each  node.  For  three  dimensions.  f.y  =  1.2,3 
associated  with  the  Jacobtans  imply  directional  identification  of  each  Jacobian  matrix  (a,,  a,,  a,.  *„  A,,  c,,. 
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L’p,  c,T,  c.,,  c,..  c.y,  Cl,,  c,..  Cl,)  With  r.  j- =  1.2,3, 4. 5  denoting  entnes  or  eacn  or  the  5  •  5  Jacooian 
matrices.  These  indices  can  be  reduced  simiiariy  tor  2-D.  Note  also  that  the  relation  i  18)  may  oe  usea  to  move 
the  delta  source  terms  to  the  RHS  of  1 33)  by  lagging  tne  source  terms  one  time  step  behind.  The  FEM  equations 
given  by  (33)  are  referred  to  as  FDMETFEM.  The  counterpan  of  (33)  based  on  FDM  schemes  of  (29)  and  (30) 
is  to  be  called  FDMETFDM. 

It  is  imponant  to  realize  that  the  integration  by  pans  as  applied  to  the  generalized  Galerkin  approximations  in 
FEM  produces  all  Neumann  boundary  integrals.  It  is  particularly  advantageous  that  Neumann  boundar>' 
conditions  through  re-evaluation  of  Jacobians  normal  to  the  boundary'  surfaces  can  simply  be  added  to  the 
boundary  nodes  for  the  stiffness  matrix  in  (35).  On  the  other  hand,  all  Neumann  boundary  conditions 
which  appear  in  (37)  act  as  source  terms.  These  features  are  absent  in  FDM.  but  implementations  of  Neumann 
boundary  conditions  can  be  handled  by  devising  special  forms  of  finite  differences  at  boundary  nodes. 

Similar  results  are  obtained  either  by  FDM  or  FEM  with  accuracy  of  computations  derived  primarily  from  the 
FDMEI  equations  of  (12).  However,  with  the  increase  of  Reynolds  number  (say  around  Re,»  10*),  it  is 
possible  that  accuracy  may  increase  with  applications  of  special  functions  such  as  Legendre  polynomials  of  high 
degree  modes  characterizing  extremely  small  turbulent  microscales.  Implementation  of  such  high  frequency 
modes  can  be  achieved  by  placing  these  modes  between  the  comer  nodes  of  isoparametric  finite  elements. 
Adaptively,  such  high  modes  can  be  chosen  as  needed  for  the  resolution  of  turbulent  microscales.  Once  again, 
the  diffusion  implicitness  parameter  5,  will  play  a  crucial  factor  in  determining  the  required  degrees  of  Legendre 
polynomial.  The  use  of  Legendre  polynomial  spectral  modes  supenmposed  onto  isoparametnc  elements  has 
been  discussed  in  [15].  Its  merit,  however,  has  not  yet  been  fully  established  for  general  applications. 

For  turbulent  flows  with  an  extremely  high  Reynolds  numoer,  the  phase  error  of  the  shon  wavelengths  can  be 
very  large.  In  this  case,  it  is  necessary  to  add  numerical  dissipation  terms  to  damp  out  the  shon  wavelengths. 
Such  numerical  viscosities  are  conceptually  different  from  the  second-order  implicitness  parameters  whose  role 
is  to  ensure  stable  solutions  while  preserving  the  solution  accuracy  dictated  by  the  first-order  implicitness 
parameters.  Toward  this  end,  it  is  desirable  to  revise  (18)  in  the  form 


At/" 


=  -2"  -Q" 


(38) 


where  Q"  is  the  numerical  dissipation  vector  in  terms  of  the  second-order  tensor  of  numerical  dissipation,  5,^. 
associated  with  the  second-order  derivatives  of 


Q"  =  K 


d^U" 
dX:  dx, 


=  fJL  Ax,  Ac 


(39) 


with  fjL  being  the  numerical  dissipation  constant  chosen  as  0  where  ^4)  is  set  approximately  equal  to  2. 

but  adjusted  from  numerical  experiments.  Note  that  the  Galerkin  integral  of  (39)  (integration  by  pans  once) 
leads  to  the  first  derivative  of  the  trial  and  test  functions  combined  wi  ij  the  nodal  values  of  U*p.  In  addition,  note 
that  the  damping  provided  by  the  second-order  derivatives  with  not  disrupt  the  formal  accuracy  of  the  FDMEI 
scheme.  This  process  may  be  applied  to  (13)  as  well. 

.  One  of  the  most  significant  aspects  of  the  FDMEI  scheme  is  that  for  low  Mach  numbers  (incompressible 
flow)  the  scheme  will  automatically  adjust  itself  to  prevent,  pressure  oscillations.  This  adjustment  is  analogous  to 
the  pressure  correction  scheme  employed  for  incompressible  flows.  Otherwise,  the  FDMEI  scheme  is  capable  of 
shock  wave  resolutions  at  high  Mach  numbers,  and  particularly  well  suited  for  dealing  with  interactions  between 
shock  waves  and  turbulent  boundary  layers  where  regions  of  high  and  low  Mach  numbers  and  Reynolds 
numbers  coexist.  In  this  case,  the  inviscid  and  viscous  interactions  are  allowed  to  take  place.  To  this  end  the 
second-order  implicitness  parameters  play  the  role  of  artificial  viscosity  needed  for  shock  wave  resolutions  in 
the  presence  of  flow  diffusion  due  to  physical  viscosity. 

In  order  to  understand  how  the  FDMEI  scheme  handies  computations  involving  both  compressible  and 
incompressible  flows,  fundamental  definitions  of  pressure  must  be  recognized.  Consider  in  the  following  that  the 
fluid  is  a  perfect  gas  and  that  the  total  energy  is  given  by 


P  1 


(40) 
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The  momentum  equation  ror  steaay  state  mcomoressinie  rotational  riow  may  oe  inteerated  to  'jive 

J  ^  3  “  J  -  -r  L‘  . )  -  .t;  a>.  j  dv 

1 

/7  -H  “  pv.v,  =  p,,  *r  w 


with 

where  is  the  component  of  a  voracity  vector.  p„  is  the  constant  of  integration,  and  m  denotes  the  spatial 
dimension. 

Combining  (.40)  and  (41)  leads  to  the  following  relationship: 

Po  =  Pic,J  +  vv,-E)-W  ^42) 

If  Pn  as  given  by  (42)  remains  a  constant,  equivalent  to  a  stagnation  pressure,  then  the  compressible  flow  as 
assumed  in  the  conservation  form  oi  the  Navier-Stokes  svstem  of  equations  has  now  been  turned  into  an 
incompressible  flow,  which  is  expected  to  occur  when  the  flow  velocity  is  sufficiently  reduced  (approximately 
.1  <  0.3  for  air).  Thus.  (42)  may  serve  as  an  equivalent  equation  of  state  for  an  incompressible  flow.  This 

can  be  identified  element  by  element  for  the  entire  domain.  Note  that  conservation  of  mass  is  achieved  for 
incompressible  flows  with  in  (42)  being  constant,  thus  keeping  the  pressure  from  oscillating. 

Once  the  Navier-Stokes  solution  via  FDMEI  is  carried  out  and  all  flow  variables  determined,  then  we 
compute  fluctuations  /'  of  any  variable  /. 


where  /  and  /  denote  the  Navier-Stokes  solution  and  its  time  average,  respectively.  This  process  may  be 
replaced  by  the  fast  fourier  transform  of  the  Navier-Stokes  solution.  Unsteady  turbulence  statistics  (turbulent 

kinetic  energy.  Reynolds  stresses,  and  various  energy  spectra)  can  be  calculated  once  the  fluctuation  quantities 
of  all  vanables  are  determined. 

Ut  us  summarize  why  the  FDMEI  .scheme  is  capable  of  handling  low  speed  and  high  speed  and  nonreacting 
and  reacting  compressible  and  incompressible  flows,  including  shock  waves  and  turbulent  flows:  ( 1)  How  is  the 
transition  from  incompressible  flow  to  compressible  flow  naturallv  and  automaticallv  accommodated  without 
using  two  separate  equations  or  two  separate  codes?  This  process  is  dictated  by  the  first-order  convection 
imp  1C itne^  parameter  j1  as  reflected  by  the  Mach  number  changes  and  the  expression  of  the  stagnation 
pressure.  (2)  How  is  the  shock  wave  captured?  As  the  Mach  number  increases  and  its  discontinuity  is  abrupt,  the 
s2  terms  associated  with  second  order  derivatives  together  with  squares  of  the  convection  Jacobian  provide 
a  equate  numencal  viscosities  through  second  order  derivatives,  similarly  as  the  Lax-Wendroff  scheme.  (3) 
How  IS  the  transition  from  laminar  to  turbulent  flows  naturally  and  automatically  accommodated?  This  process 
IS  governed  by  the  first-  and  second-order  diffusion  implicitness  parameters  (j3  and  j4)  as  calculated  from  the 
c  ^ges  o  ^  ^  eynoids  number.  The  terms  associated  with  s2  and  s4  are  responsible  for  fluctuations  of 
ve  cities,  with  the  values  of  these  implicitness  parameters  increasing  with  intensities  of  turbulence  in  coniuction 
wi  t  ®  •  usion  gradient  Jacobian  and  the  squares  of  the  diffusion  Jacobian.  This  process  allows  the 
avier  to  es  so  utions  to  contain  fluctuations  which  can  be  extracted  by  subtracting  the  time  averages  of  the 
Navier-St^es  solutions.  (4)  How  do  the  interactions  between  convection  and  diffusion  take  place?  Changes  of 
ac  num  ers  an  eynoids  numbers  as  reflected  by  both  convection  and  diffusion  implicitness  parameters 
close  to  the  wall  contribute  to  the  unsteadiness.  Away  from  the  wall,  they  contribute  to  the  transition  between 
incompressible  to  compressible  flows.  (5)  How  are  the  stiff  equations  arising  from  widely  disparate  reaction 
rates  o  a  c  emica  species  treated?  The  most  crucial  aspect  of  the  FDMEI  scheme  is  its  capability  to  identify 
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the  ratio  or  the  resiaent  time  to  the  recnon  time  as  caicuiatea  rrom  nve  dirYerent  aehniiions  or  tne  DamKohier 
numbers  between  tne  adjacent  nodal  points  and  time  steps  as  reriected  in  the  calculated  hrst-oraer  implicitness 
parameter,  j.,  and  the  second-order  implicitness  parameter,  These  parameters  provide  precise  degree  or 
computational  implicitness  at  every  nodal  point  and  every'  time  step,  contnbuting  to  the  determination  or' 
accurate  chemical  reactions. 


4.  Applications 

We  examine  here  various  example  problems:  la)  flow  over  a  flat  plate,  (b)  shock  wave  turbulent  boundarv' 
layer  interactions  on  a  compression  comer,  (c)  3D  duct  flows,  and  (d)  lid-driven  cavity  flow.  Linear 
isoparameteric  finite  elements  are  used  for  the  example  problems. 


(a)  Behavior  of  flowfield  dependent  implicitness  parameters  on  flat  plate 

First  of  all,  our  concern  is  to  test  the  behaviour  of  FDMEI  and  FDMEI-FEM.  Toward  this  objective  we 
examine  the  flow  over  a  flat  plate  investigated  earlier  by  Carter  [171  as  shown  in  Fig.  4(a),  The  initial  setting  for 
the  implicitness  parameters  are  determined  from  the  initial  conditions  of  the  flowfield  and  subsequently  updated 
after  each  time  step  until  the  steady  state  solution  is  reached. 

Corresponding  to  the  mesh  refinements  and  the  flowfields  at  steady  state  shown  in  Fig.  4(b-d),  the  contours 
ot  implicitness  parameters  and  are  given  in  Fig,  5.  It  is  seen  that  the  implicitness  parameters  themselves 
closely  resemble  the  flowfield.  There  arc  little  or  no  changes  in  Mach  numbers  and  Reynolds  numbers  between 
adjacent  nodes  or  elements  far  away  from  the  surface  of  the  plate  as  indicated  by  j  j  =  53  =  0.  Along  the  leading 
edge  shock  and  boundary  layers,  both  s^  and  move  toward  unity  indicating  that  gradients  of  all  variables 
increase.  The  final  flowfields,  as  shown  in  Fig.  4(b-d),  are  the  consequence  of  these  implicitness  parameters. 
The  implicitness  parameters  s^  and  s^  are  the  compliances  of  j,  and  53,  respectively,  with  their  primary  roles 
being  the  artificial  viscosity.  Thus,  the  first-order  implicitness  parameters  (^,,.53)  help  to  resolve  the  high 
gradients  ensuring  the  accuracy  of  the  solution.  While  on  the  other  hand,  the  second-order  implicitness 
parameters  (s^.sj  ensure  computational  stability. 

Computations  of  wall  pressure,  wall  skin  friction,  u-velocity,  t;-velocity,  density  and  temperature  distribution 
are  shown  in  Fig.  6(a-f).  The  comparison  with  the  Carter*  s  data  indicates  reasonable  agreements. 


(b)  Supersonic  flow  on  a  compression  corner 

In  this  example  we  demonstrate  calculations  of  supersonic  flow  on  a  compression  comer.  The  inlet  boundary 
conditions  inon-dimensionalized)  are  p  =  1.  M  =  2.25,  p  =  0.14,  Re  =  10^,  Pr  =  0.72,  and  i;  =  0.  with  adiabatic 
wall  condition.  The  steady  state  background  mean  flowfields  for  the  compression  comer  are  shown  in  Fig.  7(a). 
In  these  calculations,  all  perturbation  (fluctuation)  variables  are  determined  from  time  averages  of  the 
Navier-Stokes  solutions  according  to  (/^3).  The  horizontal  and  vertical  perturbation  velocities  at 

locations  close  to  the  wall  {x  =  0.10256  m,  y  =  0.001  m)  and  away  from  the  wall  (jc  =  0.10256  m,  y  =  0.04  m) 
are  shown  in  Fig.  7(b).  Note  that  u  is  extremely  unsteady  whereas  u'  is  significantly  less  unsteady  close  to  the 
wall.  Away  from  the  wall,  both  u  and  v*  are  almost  steady.  These  trends  are  reflected  in  the  turbulence 
(Reynolds)  stresses  as  shown  in  Fig.  7(c).  Turbulent  kinetic  energy  distributions  at  the  locations  upstream  of  the 
comer  {x  =  0.0513  m)  and  downstream  of  the  comer  (jc  =  0. 1333  m)  are  shown  in  Fig.  7(d).  We  observe  that  the 
turbulent  kinetic  energy  downstream  of  the  comer  is  significantly  larger  than  the  upstream.  No  turbulent 
statistics  calculations  (wave  numbers  or  frequencies  vs.  power  spectral  density)  are  attempted  at  this  time  as 
turbulence  microscales  are  not  resolved  in  this  example. 

It  should  be  noted  that  the  above  results  obtained  without  turbulence  models  or  without  the  standard  DNS 
solutions  (neither  spectral  nor  DNS-mesh  refinements)  are  regarded  as  the  consequence  of  the  time-averaging  of 
the  FDMEI  Navier-Stokes  solutions.  This  implies  that  the  fluctuation  of  variables  between  nodal  points  (Fig. 
2(a))  and  between  time  steps  (Fig.  2(b))  as  reflected  in  terms  of  the  implicitness  parameters  (5,  )  have  contributed 
to  these  physical  phenomena,  with  compressibility  and  shock  waves  dictated  by  the  Mach  number-dependent  5,, 
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'  f  problem  initial  and  adaptive  meshes  and  their  corresponding  density  contours:  (a)  Geometry  and  boundary  conditions  of 

problem.  A/,-3,  Re^  =  1000.  7*,  =  390®R;  (b)  initial  mesh  (816  elements,  875  nodes)  and  the  corresponding  density 
contours,  (max^  .  1  Kg/m  ,  min  —0.5  Kg/m  ):  (c)  one-level  adaptive  mesh  (1755  elements,  1889  nodes)  and  the  corresponding  density 

contours  (max  „  Kg/m  j;  (d)  two-level  adaptive  mesh  (4257  elements.  4547  nodes)  and  the  corresponding  density 

contours,  (max  =  2.0  Kg/m  ,  mm  =  0.5  Kg/m'). 
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Fig.  5.  Flowrield-dcpendent  Hrst-order  convection  and  diffusion  implicitness  parameters;  (a)  s^  contours:  ib)  r,  contours. 


and  with  incompressibility  and  turbulent  fluctuation  dictated  by  the  Reynolds  number  or  Peclet  number- 
dependent  (i- An  equal  participation  of  5,  and  will  be  responsible  for  shock  wave  turbulent  boundary  layer 
interactions.  A  comparison  of  the  results  of  the  FDMEI  scheme  with  the  K-e  turbulent  model  and  experimental 
data  is  shown  m  Fig.  7(e).  It  is  seen  that  the  FDMEI  results  compare  more  favorablv  with  those  of  measurements 
[18]. 

(c)  FDMEI  analysis  of  three-dimensional  flows 

To  demonstrate  the  effectiveness  of  the  flowfield-dependent  implicitness  parameters  in  3-D  flows  at  the  steady 
state,  we  examine  the  spatially  evolving  boundary  layer  (Fig.  8(a-e).  Note  that  the  contours  of  s^  and  5,  (Fig. 
8(c))  show  the  boundary  layer  effects  in  which  both  5,  and  s^  are  indicative  of  rapid  changes  of  Mach  numbers 
and  Reynolds  numbers;  respectively,  larger  (close  to  unity)  on  the  wall,  but  small  (closer  to  zero)  away  from  the 
wall.  The  velocity  vectors  and  RMS  error  distributions  versus  interactions  are  shown  in  Fig.  8(d)  and  (e), 
respectively. 

(d)  Demonstration  of  compressibility  vs.  incompressibility 

We  ask  the  question:  Can  a  single  formulation  or  computer  program  originally  designed  for  high  speed 
compressible  flows  be  applied  to  analyze  the  low  speed  incompressible  flows  The  advantage  of  FDMEI  is  to 
respond  positively  to  this  question.  To  prove  the  point,  let  us  examine  the  lid-driven  cavity  flow  at  the  steady 
state  (Fig.  9(a-f)).  Notice  that,  for  M  =  0. 1,  density  changes  occur  closer  to  the  lid,  whereas,  for  M  =  0.01, 
density  is  constant  throughout  the  domain  (Fig.  9(e)),  corresponding  to  variable  and  constant, 

respectively  (see  Eq.  (^2)).  The  equation  of  state  for  compressible  flows  is  automatically  switched  over  to 
accommodate  the  incompressible  flows.  This  advantage  is  contrary  to  the  previous  practice  such  as  the  Table 
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Fig.  8.  Three-dimensional  spatially  evolving  boundary  layer  flow:  (a)  Schematics  of  spatially  evolving  boundary  layer  flow:  (b) 
computational  geometry  and  boundary  conditions:  to  flowrtetd-dependeni  flrsi-order  convection  and  diffusion  implicitness  parameter 
contours;  (d)  velocity  vectors:  (e)  convergence  history. 


look-up  for  the  equation  of  state  for  incompressible  flow  handled  separately  through  hyperbolic  elliptic  equation 
as  derived  from  the  continuity  equation  combined  with  the  momentum  and  energy  equations.  Comparison  of  the 
results  of  FDMEI  with  those  of  the  independent  incompressible  flow  code  of  Ghia  et  al.  [19]  are  very  favorable 
as  shown  in  Fig.  9(  f)- 


Fig.  9.  Test  of  FDMEI  scheme  to  solve  incompressible  flows-lid-driven  cavity  problem:  (a)  Geometry  and  boundary  conditions:  (b)  mesh 
configuration:  (c)  streamline  contours:  (d)  vomcity  contours:  (e)  density  distribution:  (f)  u-  and  u-vclocity  compansons  with  Ghia  ci  al.  [19]: 
(g)  stagnation  pressure  at  M  =0.0!  and  M  =0.1. 
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5.  Concluding  remarks 

The  validity  of  the  proposed  new  approach  to  computational  fluid  dynamics  has  been  demonstrated  through 
some  example  problems.  Excluded  from  these  examples  are  reacting  flows  which  are  reponed  elsewhere  [16]. 
Also  excluded  is  the  effect  of  additional  spectral  modes  of  Legendre  polynomials  which  are  described  in  [15]. 

one  of  the  example  problems  have  been  earned  out  with  mesh  refinements  required  for  resolving  turbulent 
microscales  due  to  the  limitation  of  computer  time.  The  following  concluding  remarks  are  provided: 

(a)  The  flowfield-dependent  implicitness  parameters  as  calculated  from  the  current  flowfield  information  are 
indicative  of  the  magnitude  of  gradients  of  all  variables  and  adjust  the  computational  schemes  accordingly 

or  every  nodal  point  or  element,  rather  than  dictated  by  arbitrarily  selected  constant  parameters 
throughout  the  domain. 

(b)  ^e  first-order  implicitness  parameters  j,,  Sy  and  as  calculated  from  Mach  numbers,  Reynolds  or 

eclet  numbers,  and  Damkohler  numbers,  respectively,  ensure  the  solution  accuracy,  whereas  the 
second-order  implicitness  parameters  and  5^  which  are  determined  as  compliances  of  j,,  and  x., 
respectively,  assist  in  the  solution  stability. 

(c)  The  FDMEI  method  is  capable  of  resolving  mutual  interactions  and  transition  between  viscous  and 
inviscid  flows,  compressible  and  incompressible  flows,  and  laminar  and  turbulent  flows,  in  ail  speed 
re^mes.  Further  research  on  FDMEI  is  required  in  order  to  investigate  many  other  physical  phenomena 
including  hypersonic  and  high  temperature  flows  in  3D. 


Appendix  A.  Analogies  of  FDMEI  to  currently  available  FDM  and  FEM  schemes 
Analogies  of  FDMEI  to  currently  available  computational  schemes  of  FDM  and  FEM  are  summarized  below. 

A./.  Analogies  of  FDMEI  to  FDM 

Some  of  the  FDM  schemes  are  compared  with  FDMEI  in  Table  A.l. 
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Other  schemes  of  FDM  are  compared  with.  FDMEI  as  follows: 

(a)  Lax-Wendrojf  scheme 


The  Lax-Wendroff  scheme  without  artificial  viscosity  takes  the  form 

This  scheme  arises  if  we  set  in  FDMEI 

5j”0,  S2  ^  0  ,  0,  0 

(b)  Lax-Wendrojf  scheme  with  viscosity 
The  Lax-Wendroff  scheme  with  viscosity  is  given  by 


(A.1) 


(A.2) 


with 


At 

^/  +  i/2  2  2  Ax  ^ +  i/2^^/  +  i 

♦  Ff  +  F,^.  At 

This  scheme  arises  if  we  set 


I  /  2  -  I  /  2  !♦  5'2~0»  Sy  —  0,  54  —  0 

This  implies  that  aj,  in  FDMEI  plays  a  role  of  artificial  viscosity  which  is  manually  implemented  in  the 
Lax-Wendroff  scheme. 


(c)  Explicit  McCormack  scheme 

Combining  the  predictor  corrector  steps  of  McCormack  scheme  we  write 
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The  FDMEI  becomes  identical  to  this  scheme  with  the  toiiowing  adjustments: 


(A.3) 


^/+l/2  ^-1/2  ^ 

5,  =  0.  •^2=0.  .^3  =  0,  =0 

and  the  j.  term  in  the  FDMEI  method  is  equivalent  to 

^/=f  ^6u:-wi, -u%,) 

This  again  is  a  manifestation  that  shows  the  equivalent  of  the  s.  terms  is  manually  supplied  in  the  McCormack 
method. 


(d)  First-order  upwind  scheme 


This  scheme  is  written  as 
A/ 


=  - {  [t  ^  f/;* )]  -  [^  T  ^ ^ 

The  FDMEI  analogy  is  obtained  by  setting 
1 


F"  =  —  F" 

^  i  2  ^  ' 


F"  =  — F" 

' "  *  2  ^  • 

-2Af/::;  ^A(;:*:,')  =  ia|{£/;\, 

where  C  is  the  Courant  number. 


(e)  Implicit  McCormack  scheme 

With  all  second-order  derivatives  removed  from  (11)  we  obtain  the  implicit  McCormack  Scheme  by  setting 
s^  =  1,  Sy=0,  s^=0.  However,  it  is  necessary  to  divide  the  process  into  the  predictor  and  corrector 

steps.  Once  again,  the  flowfieid-dependent  implicitness  parameters  for  FDMEI  will  allow  the  computation  to  be 
performed  in  a  single  step. 

(f)  PISO  and  SIMPLE 

The  basic  idea  of  PISO  and  SIMPLE  is  analogous  to  FDMEI-FEM  in  that  the  pressure  correction  process  is  a 
separate  step  in  PISO  or  SIMPLE,  whereas  the  concept  of  pressure  correction  is  implicitly  embedded  in 
FDMEI-FEM  by  updating  the  implicitness  parameters  based  on  the  upstream  and  downstream  Mach  numbers 
and  Reynolds  numbers  within  an  element. 

The  elliptic  nature  of  the  pressure  Poisson  equation  in  the  pressure  correction  process  resembles  the  terms 
embedded  in  the  terms  in  (3£-).  Specifically,  examine  the  terms  involving  and  and 

term  involving  All  of  these  terms  are  multiplied  by  which  provide  dissipation  against  any 

pressure  oscillations.  Question:  Exactly  when  is  such  dissipation  action  needed?  This  is  where  the  importance  of 
implicitness  parameters  based  on  flowfield  parameters  comes  in.  As  the  Mach  number  becomes  very  small 
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I incomoressibility  erfects  dominate)  the  impiicitness  parameters  5,  «nd  s,  jaiculated  from  tne  current  riowneia 
will  be  indicative  of  pressure  correction  required.  Notice  mat  a  aeiicate  oaiance  oeiween  Mach  numoer  [s.  i> 
Mach  number  dependent)  and  Reynolds  numoer  or  Peclei  number  ia,  Reynolds  number  or  Peciet  number 
dependent)  is  a  crucial  factor  in  achieving  convergent  and  stable  solutions.  Of  course,  on  the  other  hand,  high 
Mach  number  flows  are  also  dependent  on  these  implicitness  parameters.  In  this  case  all  implicitness  parameters, 
i,.  5.,  .y,,  will  play  important  roles. 


A.2.  Analogies  of  FDMEl  to  FEM 

(a)  Generalized  Taylor  Galerkin  iGTG)  with  convection  and  diffusion  Jacobians 

Earlier  developments  for  the  solution  of  Navier-Slokes  system  of  equations  were  based  on  GTG  without 
using  the  implicitness  parameters.  They  can  be  shown  to  be  special  cases  of  FDMEI-FEM. 

In  terms  of  both  diffusion  Jacobian  and  diffusion  gradient  Jacobian,  we  write 

ac,  ^  hu 


hG  flG  hU 

^=“77-  V^=  — 
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Thus,  It  follows  from  ( 10) 


with  5,  =  s^  =  .Tj  =  s^  =  s^  =  Q  and  5,  “  ^ 
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Using  the  definitions  of  convection,  diffusion,  and  diffusion  rate  Jacobians  discussed  in  Section  2,  the  temporal 
rates  of  change  of  convection  and  diffusion  variables  may  be  written  as  follows: 
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Substituting  (A.6)  and  (A.7)  into  (A.5)  yields 
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and  neglecting  the  spatial  and  temporal  derivatives  of  we  rewrite  (A. 8)  in  the  form 
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Here,  the  .second  derivatives  of  G,  are  neglected  and  all  Jacoo.ans  are  assumea  to  remain  constant  w.thin  an 
incremenial  time  step  but  updated  at  subsequent  time  steps. 

Applying  the  Galerkin  rinite  element  formulation,  we  have  an  implicit  scheme. 

ag;:  '  =  h:^  ^  .v;:: '  - 
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Here.^we  note  that  the  algorithm  given  by  (A.IO)  results  from  (i3)  bv  seitine  J,  =J,  =5,=0  s  =1 
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(b)  CrC  inV/i  convection  Jacobians 

Diffusion  Jacobians  may  be  neglected  if  their  influences  is  negligible.  In  this  case  the  Taylor-Galerkin  rtnite 
.tem=„,  .„a;og  may  ba  darned  using  only  ,be  convecve  Jacob, an  from  the  Taylor  senes  expanstT 
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Substituting  (A.12)  and  (A.13)  into  (A.ll).  we  obtain 
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where  second  derivatives  of  G.  is  assumed  to  be  negligible  and  B  is  constant  in  space  and  time, 
implicit  finite  element  scheme. 
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It  should  be  noted  that  the  form  (A.1 5  )  arises  from  (23)  with  5,  -  =  s^-  b^-Q  and  j.  = 

similar  to  Hassan  et  al.  [13]. 

(c)  Generalized  Petrov -Giilerkin  {GPG) 

The  Generalized  Petrov-Galerkin  (GPG)  method  can  be  identified  by  setting  =5.= 
.=  rf  =  0.  Q”  =  0.  and  =  3  so  that  ( 10)  takes  the  form 
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For  the  steady  state  non-incremental  form  in  l-D  we  write  (A.  16)  in  the  form 
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Taking  the  Galerkin  integral  of  (A.  17)  leads  to 
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with  a  =  C/2  and  C  =  i;  St/Xx  being  the  Courant  number. 

For  isoparametric  coordmate.s  m  two  uimensions.  the  Petrov-Galerkin  test  tunction  assumes  me  torm 

(A.20) 
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(d)  Cl.araaenmc-hased  Ziankwmcz-Codincr  scheme 
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Continuity 

The  continuity  equation  can  be  extracted  from  (A.21)  by  setting  as  follows: 
At/""'  ^Ap"*' 

F"  — ^  pv'  — >  U" 

s^a^  AU"  '  —  V|  Ap  u"  ~ri,  AU" 

i/2a,F"  6,^ 

I  /2y,a,«  At/"" '  61,  ly,  Ap"  ■ ' 

These  substitutions  to  (A.2I)  lead  to 

Ap""'=(jTAp)''"  =  -Arf’'’^^"'>'' 


(A.21) 


d(ApV^)" 


\dx,d.x,  dx,dx^ 


(A.22) 
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which  is  identical  to  Eq.  (33)  in  (141  with  (Am,  i"*'  u  •  ..  ■  '  -J 

correction  process  in  case  the  flow  is  incompressible.  'n'ermediate  step.  This  represents  the  pressure 

Momentum 

A  similar  procedure  can  be  applied  to  (A  2n  for  rho 

yyiKu  lo  ror  the  momentum  equations. 
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which  is  equivalent  to  Eq.  (30)  in  (14]  with  a  =  f  .  1  -  =  5,,  and  all  terms  of  5^,,  and  c  ,  being  negiectea  m 

FDMEL 


Energy 

Again,  from  (A,21),  neglecting  ail  n  -r  1  terms  of  the  FDMEI  equauon.  we  obtain 
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d{pEv.  +  vpY 
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(A.24) 


which  corresponds  to  Eq.  (40)  of  [14].  The  solution  steps  begin  with  (A.23),  followed  by  (A.21)  and  (A.24),  and 
continue  iteratively  until  convergence.  Note  that  the  pressure  corrections  for  incompressibility  are  internally 
carried  out  in  FDMEI  as  the  pressure  second  derivatives  arise  in  Eq.  (10).  Note  also  that  in  FDMEI  all  implicit 
terms  may  be  retained  for  computational  accuracy  and  efficiency  for  any  physical  situation. 
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